---
output: 
  bookdown::pdf_document2:
    toc: false
    citation_package: natbib
    fig_caption: yes
    latex_engine: pdflatex
    template: ~/Dropbox/Harvard/pandoc-default-latex.tex
    keep_tex: false
always_allow_html: true
title: "Replication Package: The Politics of Disaster Prevention"
author: 
- name: "Martin Gilens"
  affiliation: "University of California, Los Angeles"
  email: gilens@ucla.edu
  orcid_id: 0000-0001-5485-6275
- name: "Tali Mendelberg"
  affiliation: "Princeton University"
  email: talim@princeton.edu
  orcid_id: 0000-0002-4494-7541
- name: "Nicholas Short"
  affiliation: "Princeton University"
  email: nick.short@princeton.edu
  orcid_id: 0000-0002-2401-8315
date: "`r format(Sys.time(), '%B %d, %Y')`"
header-includes:
    \usepackage{mathtools}
    \usepackage{natbib}
    \usepackage{pdflscape}
    \usepackage{longtable}
    \usepackage{setspace}
    \usepackage{pdfpages}
    \usepackage{caption}
    \captionsetup[table]{font=large}
    \pagenumbering{gobble}
geometry: margin=1in
fontfamily: mathpazo
fontsize: 12pt
linestretch: 2
bibliography: manuscript.bib
biblio-style: apsr2006.bst
---
\thispagestyle{empty}

\clearpage
\pagenumbering{arabic} 

\pagebreak

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(tidyverse)
library(estimatr)
library(survey)
library(kableExtra)
library(stargazer)

options(dplyr.summarise.inform = FALSE)
```

# Notes

The code in this file, and in the accompanying Stata .do file, generates all of the figures, tables, and statistics published or discussed in the manuscript and in the online appendix for The Politics of Disaster Prevention.

When you compile (or "knit") this RMarkdown document, it will generate a PDF file of the same name (`the-politics-of-disaster.pdf`).  Currently, it does not store hard copies of the figures, but the user or replication analyst can accomplish this by changing the `keep_tex` option in the YAML header from FALSE to TRUE.  This will not only preserve PDFs of each figure, it will also generate and save a LaTex copy of the output (`the-politics-of-disaster.tex`).

The replication package contains an empty folder called `tables`.  Upon compiling, this document will generate each analytical table in the paper, save the .tex copies to that folder, and input the tables into `the-politics-of-disaster.pdf`.  That way the replication analyst can refer to the raw .tex files, but also compare the rendered PDF image to the PDFs for the manuscript and online appendix.  (Non-analytical, manually generated, tables are not reproduced herein).

The `stargazer` library for publishing nicely formatted latex tables unfortunately does not accept linear models generated with `estimatr`.  The code below works around this issue by estimating standard linear models with `lm`, extracting robust standard errors with `sandwich`, and then passing both objects to `stargazer`.  However, we use `estimatr` to quickly generate the same output for plotting regression output in figures.

\pagebreak

# Manuscript

Load the data files.

```{r, message = FALSE, warning = FALSE}

# Load the 2021 survey data

load(file = "2021-survey-data.RData")
survey_1_design <- svydesign(id = ~1, weights = ~weights, data = survey_1, na.rm = T)

# Write the 2021 survey data file to .dta for executing certain analyses in Stata

survey_1 %>% haven::write_dta("2021-survey-stata-dataset.dta")

# Load labeled data output from the textual analysis.  This contains "mentions" for all observations and "support" for the hand-coded data.

load(file = "2021-open-text-1-data.RData")
opentext_design <- svydesign(id = ~1, weights = ~weights, data = opentext_responses, na.rm = T)

# Subset the open-text design object for those who mention prevention and relief

opentext_design_mention_prev <- subset(opentext_design, mention_prev == 1)
opentext_design_mention_relief <- subset(opentext_design, mention_relief == 1)

# Load the cleaned open text responses (all lower case, punctuation removed)

load(file = "2021-open-text-1-responses.RData")
load(file = "2021-open-text-2-responses.RData")

# Load the combined survey (2021 and 2023)

load(file = "combined-survey-data.RData")

# Create a survey 2 dataset

survey_2 <- combined_survey %>% filter(survey == "survey_2")
survey_2_design <- svydesign(id = ~1, weights = ~weights, data = survey_2, na.rm = T)

# Create subsets of the data needed for regression analysis

survey_2_nd_and_hd <- survey_2 %>% filter(condition == "natural_and_health_disasters")

survey_2_nd_or_hd <- survey_2 %>%
  filter(condition %in% c("natural_disasters", "health_disasters")) %>%
  mutate(condition = factor(condition, levels = c("natural_disasters", "health_disasters")))

combined_survey_same_wording <- combined_survey %>%
  filter(is.na(condition) | condition == "natural_and_health_disasters")

combined_surveys_same_wording_design <- svydesign(id = ~1, weights = ~weights, data = combined_survey_same_wording, na.rm = T)

# Create a dataset with binary outcomes for logistic regression (testing mentions on policy preferences)

opentext_responses_binary <- opentext_responses %>%
  mutate(across(.cols = c(more_relief, more_prevention),
                .fns = ~ case_when(.x %in% c(0.75,1) ~ 1,
                                   is.na(.x) ~ NA_real_,
                                   TRUE ~ 0)),
         across(.cols = c(more_relief, more_prevention, prefer_prevention),
                .fns = factor))
```


## Figures

The code below generates Figure 1 in the manuscript.

```{r overall-support}

# Calculate mean support and standard errors for each spending category, plus 83 percent confidence intervals (lower and upper bound)

overall_support <- as_tibble(svymean(~health+highways_and_bridges+natural_and_health_disasters+more_relief+more_prevention, survey_1_design, na.rm = TRUE), rownames = "variable") %>%
  mutate(variable = case_when(variable == "health" ~ "1: Health",
                              variable == "highways_and_bridges" ~ "2: Highways and \n bridges",
                              variable == "natural_and_health_disasters" ~ "3: Natural and public \n health disasters",
                              variable == "more_prevention" ~ "4: Disaster \n prevention",
                              variable == "more_relief" ~ "5: Disaster \n relief",
                              TRUE ~ NA_character_),
         lower = mean - 1.386*SE, 
         upper = mean + 1.386*SE) 

# Set x-axis breaks and labels for plotting

x_breaks <- c(0, 0.25, 0.5, 0.75, 1.0)
x_labels <- c("Much less", "Somewhat less", "Same", "Somewhat more", "Much more")

# Plot the results with a dashed horizontal line separating broad spending categories from disaster sub-categories

overall_support %>%
  ggplot(aes(x = mean, y = reorder(variable, desc(variable)), xmin = lower, xmax = upper)) + 
  geom_point(shape = 1) +
  geom_errorbarh(aes(height = .01)) + 
  geom_hline(yintercept = 2.5, linetype = 2) +
  labs(x = NULL, y = NULL) + 
  theme_minimal() +
  scale_x_continuous(breaks = x_breaks, labels = x_labels, limits = c(0,1)) + 
  theme(axis.text.x = element_text(angle = 45, size = 12, hjust = 1), axis.text.y = element_text(size = 12))
```

\pagebreak

The code below generates Figure 2 in the manuscript.

```{r marginal-means}

# Make the outcomes dichotomous (prefer_prevention and prefer_pandemic_prevention are already binary)

marginal_means <- survey_1 %>%
  mutate(across(.cols = c(more_relief, more_prevention, more_pandemic_relief, more_pandemic_prevention, more_covid_preparation),
                .fns = ~ case_when(.x >= 0.75 ~ 1,
                                   is.na(.x) ~ NA_real_,
                                   TRUE ~ 0)))

# Create a survey design object from the re-coded data

marginal_means_svy <- svydesign(id = ~1, weights = ~weights, data = marginal_means, na.rm = T)

# Build a data set with means and SE for each subset of the survey, label sub-groups, and recode dependent variables

mm_with_se <- as_tibble(svymean(~more_relief+more_prevention+prefer_prevention+more_pandemic_relief+more_pandemic_prevention+prefer_pandemic_prevention+more_covid_preparation, 
        design = subset(marginal_means_svy, party == "Republican"), na.rm = TRUE), rownames = "dep_var") %>% mutate(sub_group = "Republican") %>%
  bind_rows(
    as_tibble(svymean(~more_relief+more_prevention+prefer_prevention+more_pandemic_relief+more_pandemic_prevention+prefer_pandemic_prevention+more_covid_preparation, 
        design = subset(marginal_means_svy, biden == "Voted for Trump"), na.rm = TRUE), rownames = "dep_var") %>% mutate(sub_group = "Voted for Trump")) %>%
  bind_rows(
    as_tibble(svymean(~more_relief+more_prevention+prefer_prevention+more_pandemic_relief+more_pandemic_prevention+prefer_pandemic_prevention+more_covid_preparation, 
        design = subset(marginal_means_svy, trust_in_gov == "Low"), na.rm = TRUE), rownames = "dep_var") %>% mutate(sub_group = "Low trust in gov.")) %>%
  bind_rows(
    as_tibble(svymean(~more_relief+more_prevention+prefer_prevention+more_pandemic_relief+more_pandemic_prevention+prefer_pandemic_prevention+more_covid_preparation, 
        design = subset(marginal_means_svy, need_experts == "Low"), na.rm = TRUE), rownames = "dep_var") %>% mutate(sub_group = "Low need for experts"))  %>%
  mutate(dep_var = case_when(dep_var == "more_relief" ~ "Spend more on \nrelief",
                             dep_var == "more_prevention" ~ "Spend more on \nprevention",
                             dep_var == "prefer_prevention" ~ "Prioritize \nprevention",
                             dep_var == "more_pandemic_relief" ~ "Spend more \npand. relief",
                             dep_var == "more_pandemic_prevention" ~ "Spend more \npand. prev.",
                             dep_var == "prefer_pandemic_prevention" ~ "Prioritize \npand. prev.",
                             dep_var == "more_covid_preparation" ~ "Spend more \nCovid prep.",
                             TRUE ~ NA_character_))

# Calculate the number of observations in each sub-group and create annotations using these n's

rep <- marginal_means %>%
  filter(party == "Republican") %>%
  summarize(n = n()) %>%
  mutate(sub_group = "Republican")

trump <- marginal_means %>%
  filter(biden == "Voted for Trump") %>%
  summarize(n = n()) %>%
  mutate(sub_group = "Voted for Trump")

low_trust_in_gov <- marginal_means %>%
  filter(trust_in_gov == "Low") %>%
  summarize(n = n()) %>%
  mutate(sub_group = "Low trust in gov.")

low_need_experts <- marginal_means %>%
  filter(need_experts == "Low") %>%
  summarize(n = n()) %>%
  mutate(sub_group = "Low need for experts")

mm_n_obs <- rep %>%
  bind_rows(trump) %>%
  bind_rows(low_trust_in_gov) %>%
  bind_rows(low_need_experts) 

annotations <- mm_n_obs %>%
  select(sub_group, n) %>%
  distinct() %>%
  mutate(x = rep(0.75, 4),
         y = c(1:4),
         label = paste("n = ", as.character(n), sep ="")) %>%
  as.data.frame()

# Figure 2

mm_with_se %>%
  filter(dep_var %in% c("Spend more on \nprevention","Spend more on \nrelief","Prioritize \nprevention")) %>%
  mutate(dep_var = factor(dep_var),
         x_min = mean - 1.386*SE,
         x_max = mean + 1.386*SE) %>%
  ggplot(aes(x = mean, y = sub_group, fill = dep_var, xmin = x_min, xmax = x_max)) + 
  geom_bar(width = 0.5, stat = "identity", position = position_dodge(c(0.5, 0.5, 0.7))) +
  geom_errorbar(width = 0.1, position = position_dodge(c(0.5, 0.5, 0.7))) + 
  theme_minimal() + 
  labs(x = NULL, y = NULL, fill = "Outcome", size = 11) + 
  xlim(c(0,1)) + 
  geom_vline(xintercept = 0.5, linetype = "dashed") + 
  geom_text(data = annotations, aes(x = x, y = y, label = label, fill = NULL, xmin = NULL, xmax = NULL), size = 4) + 
  theme(legend.position = "bottom", 
        axis.text.y = element_text(size = 11), 
        legend.text = element_text(size = 11)) +
  scale_fill_manual(values = c("#00BA38","#619CFF","#F8766D"), breaks = c("Spend more on \nrelief", "Spend more on \nprevention", "Prioritize \nprevention"))
```

\pagebreak

The code below generates Figure 3 in the manuscript.

```{r reasons}

# Calculate mean persuasiveness for all prevention reasons

prev_reasons <- as_tibble(svymean(~Q8_1+Q8_2+Q8_3+Q8_4+Q8_5+Q8_6+Q8_7+Q8_8+Q8_9, survey_1_design, na.rm = TRUE), rownames = "variable") %>%
  mutate(lower = mean - 1.386*SE,
         upper = mean + 1.386*SE,
         group = "Reasons to prioritize prevention")

# Calclulate mean persuasiveness for all relief reasons

relief_reasons <- as_tibble(svymean(~Q11_1+Q11_2+Q11_3+Q11_4+Q11_5+Q11_6+Q11_7+Q11_8+Q11_9+Q11_10, survey_1_design, na.rm = TRUE), rownames = "variable") %>%
  mutate(lower = mean - 1.386*SE,
         upper = mean + 1.386*SE,
         group = "Reasons to prioritize relief")

# This function will convert the variable numbers to short synopses of each reason text

variable_cleaning <- function(x){
  x <- case_when(x == "Q8_1" ~ "D: Experts say disasters very likely to happen again",
                 x == "Q8_2" ~ "A: Disasters cause enormous personal losses",
                 x == "Q8_3" ~ "F: Disasters can result in drastic measures",
                 x == "Q8_4" ~ "G: Not enough time to save many victims after disaster happens",
                 x == "Q8_5" ~ "B: Some experts say $1 on prevention saves $15 in losses",
                 x == "Q8_6" ~ "H: Steady investment is better than sporadic relief",
                 x == "Q8_7" ~ "C: Prevents suffering by the most vulnerable",
                 x == "Q8_8" ~ "E: Avoids harm to front-line workers (black, white, Hispanic, and others)",
                 x == "Q8_9" ~ "I: A lot of benefit will go to immigrants and minorities",
                 x == "Q11_1" ~ "J: Some experts say severe disasters not likely to happen again soon",
                 x == "Q11_2" ~ "A: Disasters cause enormous personal losses",
                 x == "Q11_3" ~ "G: There is enough time to save many victims after disaster happens",
                 x == "Q11_4" ~ "D: Relief spending is transparent, prevention is not",
                 x == "Q11_5" ~ "B: Without relief, you lose ability to send help to the right place",
                 x == "Q11_6" ~ "F: Relief goes to identifiable individuals",
                 x == "Q11_7" ~ "E: You know who is getting help and if they are deserving",
                 x == "Q11_8" ~ "H: Better to meet urgent current needs",
                 x == "Q11_9" ~ "C: Can't prevent all disasters and we need money for relief",
                 x == "Q11_10" ~ "I: A lot of money will go to immigrants and minorities",
                 TRUE ~ NA_character_
                 )
}


# Join the datasets together

reasons <- prev_reasons %>%
  bind_rows(relief_reasons) %>%
  mutate(group = factor(group))

# Figure 3

x_breaks <- c(0, 0.33, 0.66, 1.0)
x_labels <- c("Not a good \n reason", "A fair \n reason", "A good \n reason", "A very good \n reason")
reasons %>%
  mutate(variable = variable_cleaning(variable)) %>%
  ggplot(aes(x = mean, y = reorder(variable, desc(variable)))) +
  geom_point(shape = 1, size = 1) +
  geom_errorbarh(aes(xmin = lower, xmax = upper, height = .1)) +
  theme_bw() +
  labs(x = NULL, y = NULL) +
  facet_wrap(~group, nrow = 2, ncol = 1, scale = "free_y") +
  scale_x_continuous(breaks = x_breaks, labels = x_labels, limits = c(0,1)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        axis.text.y = element_text(size = 10),
        legend.text = element_text(size = 11))
```

\pagebreak

Note, Figure 4 in the manuscript is a flowchart that was created manually.

\pagebreak

The code below generates Figure 5 in the manuscript.

```{r stability}

# Subset the survey design object to include only those who saw prevention reasons first (and relief second) and only those who saw relief reasons first (and prevention second)

prev_first <- subset(survey_1_design, reasons_prev_first == 1)
relief_first <- subset(survey_1_design, reasons_prev_first == 0)

# Calculate average levels of support (on five point scale) in both subsets and join the datasets together

avg_priorities <- as_tibble(svymean(~after_reasons_prefer_prevention+after_reasons_prefer_relief, prev_first, na.rm = TRUE), rownames = "variable") %>%
  mutate(
    group = case_when(variable == "after_reasons_prefer_prevention" ~ "Reasons from One Side",
                      variable == "after_reasons_prefer_relief" ~ "Reasons from Both Sides")
    ) %>%
  bind_rows(
    as_tibble(svymean(~after_reasons_prefer_relief+after_reasons_prefer_prevention, relief_first, na.rm = TRUE), rownames = "variable") %>%
      mutate(
        group = case_when(variable == "after_reasons_prefer_relief" ~ "Reasons from One Side",
                          variable == "after_reasons_prefer_prevention" ~ "Reasons from Both Sides")
        )) %>%
  mutate(
    lower = mean - 1.386*SE,
    upper = mean + 1.386*SE
    ) %>%
  mutate(group = factor(group, ordered = TRUE, levels = c("Reasons from One Side", "Reasons from Both Sides"))) %>%
  mutate(variable = case_when(variable == "after_reasons_prefer_prevention" ~ "Prioritize prevention over relief",
                              variable == "after_reasons_prefer_relief" ~ "Prioritize relief over prevention"))


# Figure 5

x_breaks <- c(0, 0.25, 0.5, 0.75, 1.0)
x_labels <- c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Strongly agree")

avg_priorities %>%
  ggplot(aes(x = mean, y = reorder(variable, desc(variable)), xmin = lower, xmax = upper)) +
  geom_point(shape = 1, size = 1.2) + # aes(shape = group), position = position_dodge(width = 0.5)
  geom_errorbarh(height = 0.02) + # position = position_dodge(width = 0.5)
  theme_bw() +
  scale_x_continuous(breaks = x_breaks, labels = x_labels, limits = c(0,1)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = NULL, y = NULL) +
  theme(legend.position = "bottom",
        axis.text.y = element_text(size = 11),
        axis.text.x = element_text(size = 11)) +
  facet_wrap(~group, scales = "free_x")
```

\pagebreak

The code below produces long-format (multiple observations per respondent) datasets for doing regression analysis on the candidate evaluations.

```{r cand-eval-long-datasets}

# Clean datasets for analyzing candidate choice experiment (pivot to long format, with each respondent evaluation four candidates)

candidate_data <- survey_1 %>%
  filter(!is.na(cand_voted)) %>% # Remove those who were not assigned to one of the four conditions
  select(ResponseId, cand_more_prev, cand_not_prev, cand_more_relief, cand_not_relief, cand_voted, weights) %>%
  pivot_longer(-c(ResponseId, cand_voted, weights), names_to = "condition", values_to = "support") %>%
  filter(!is.na(support)) %>%
  mutate(condition = factor(condition, levels = c("cand_not_relief","cand_not_prev","cand_more_relief","cand_more_prev")))

# Construct a similar dataset, but limit to first evaluations

candidate_first_data <- survey_1 %>%
  filter(!is.na(cand_voted), !is.na(first_cand_more_prev) | !is.na(first_cand_more_relief)) %>% # Remove those who were not assigned to one of the four conditions and non-first candidates
  select(ResponseId, first_cand_more_prev, first_cand_more_relief, cand_voted, weights) %>%
  pivot_longer(-c(ResponseId, cand_voted, weights), names_to = "spendingtype", values_to = "support") %>%
  filter(!is.na(support)) %>%
  mutate(spendingtype = factor(spendingtype, levels = c("first_cand_more_relief","first_cand_more_prev")))

# For plotting reverse the contrasts to easily extract a single coefficient

# Set the contrasts so that relief spending / did not increase spending is the baseline condition

candidate_relief_contrast <- candidate_data %>%
  mutate(spendingtype = ifelse(condition %in% c("cand_more_prev","cand_not_prev"), "prevention", "relief"),
         spendingposition = ifelse(condition %in% c("cand_more_prev","cand_more_relief"), "increased", "did not increase"),
         spendingtype = factor(spendingtype, levels = c("relief","prevention")),
         spendingposition = factor(spendingposition, levels = c("did not increase", "increased")))

# Set the contrasts so that prevention spending / did not increase spending is the baseline condition

candidate_prev_contrast <- candidate_data %>%
  mutate(spendingtype = ifelse(condition %in% c("cand_more_prev","cand_not_prev"), "prevention", "relief"),
         spendingposition = ifelse(condition %in% c("cand_more_prev","cand_more_relief"), "increased", "did not increase"),
         spendingtype = factor(spendingtype, levels = c("prevention","relief")),
         spendingposition = factor(spendingposition, levels = c("did not increase","increased")))

# Do the same for the first-candidate-only dataset

candidate_first_prev_contrast <- candidate_first_data %>%
  mutate(spendingtype = case_when(spendingtype == "first_cand_more_prev" ~ "prevention",
                                  spendingtype == "first_cand_more_relief" ~ "relief",
                                  TRUE ~ NA_character_),
         spendingtype = factor(spendingtype, levels = c("prevention","relief"))
         )

candidate_first_relief_contrast <- candidate_first_data %>%
  mutate(spendingtype = case_when(spendingtype == "first_cand_more_prev" ~ "prevention",
                                  spendingtype == "first_cand_more_relief" ~ "relief",
                                  TRUE ~ NA_character_),
         spendingtype = factor(spendingtype, levels = c("relief","prevention"))
         )
```

The code below generates Figure 6 in the manuscript.

```{r candidate-selection-barplot}

# Run regressions again (using lm_robust for plotting)

candidate_relief_results <- lm_robust(support ~ cand_voted + spendingposition * spendingtype, 
                                      data = candidate_relief_contrast, 
                                      weights = weights, 
                                      clusters = ResponseId, 
                                      se_type = "stata")

candidate_prev_results <- lm_robust(support ~ cand_voted + spendingposition * spendingtype, 
                                    data = candidate_prev_contrast, 
                                    weights = weights, 
                                    clusters = ResponseId, 
                                    se_type = "stata")

# Keep desired estimates and bind together

relief_did_figure <- candidate_relief_results %>%
  tidy() %>%
  filter(term == "spendingpositionincreased") %>%
  mutate(term = ifelse(term == "spendingpositionincreased", "Increased spending", term),
         model = "relief")

prev_did_figure <- candidate_prev_results %>%
  tidy() %>%
  filter(term == "spendingpositionincreased") %>%
  mutate(term = ifelse(term == "spendingpositionincreased", "Increased spending", term),
         model = "prevention")

did_figure <- relief_did_figure %>%
  bind_rows(prev_did_figure)

# Figure 6

did_figure %>%
  mutate(model = case_when(model == "relief" ~ "Relief",
                           model == "prevention" ~ "Prevention"),
    term = "") %>% # Doing this drops what looks like a y-axis label that says "Increased spending"
  ggplot(aes(x = estimate, y = model, fill = model)) + 
  geom_bar(width = 0.2, stat = "identity", position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(xmin = estimate-1.96*std.error, xmax = estimate+1.96*std.error), width = 0.05, position = position_dodge(width = 0.5)) + 
  theme_minimal() + 
  theme(legend.position = "none",
        axis.text = element_text(size = 12),
        axis.title.x = element_text(size = 12)) + 
  labs(x = "Estimated increase in likelihood of voting for candidate", y = NULL) + 
  scale_fill_manual(values = c("#00BA38","#619CFF"), breaks = c("Relief", "Prevention"))
```

\pagebreak

The code below generates Figure 7 in the manuscript.

```{r open-text-bar}

# Calculate the share of all respondents (n = 2,932) who mentioned each concept; in the hand-coded sample (n = 584), in the subset of respondents who mentioned each concept, calculate the share who expressed support.  Bind the estimates together, clean the variable names, create a binary variable for spending type (relief or prevention), and calculate 83 percent confidence intervals from standard errors.

opentext_means <- as_tibble(svymean(~mention_relief+mention_prev, opentext_design, na.rm = TRUE), rownames = "variable") %>%
  bind_rows(as_tibble(svymean(~pos_relief, opentext_design_mention_relief, na.rm = TRUE), rownames = "variable") %>% rename(SE = pos_relief)) %>%
  bind_rows(as_tibble(svymean(~pos_prev, opentext_design_mention_prev, na.rm = TRUE), rownames = "variable") %>% rename(SE = pos_prev)) %>%
  mutate(variable = case_when(variable == "mention_relief" ~ "Mentions relief",
                              variable == "mention_prev" ~ "Mentions prevention",
                              variable == "pos_relief" ~ "Support for relief",
                              variable == "pos_prev" ~ "Support for prevention"),
         spending_type = case_when(variable == "Mentions relief" ~ "Relief",
                                   variable == "Mentions prevention" ~ "Prevention",
                                   variable == "Support for relief" ~ "Relief",
                                   variable == "Support for prevention" ~ "Prevention"),
         lower = mean - 1.386*SE, 
         upper = mean + 1.386*SE,
         variable = factor(variable, levels = c("Support for prevention", "Support for relief", "Mentions prevention", "Mentions relief")),
         spending_type = factor(spending_type, levels = c("Relief","Prevention"))
         )

# Plot the shares and confidence intervals, using spending_type to color the bars

opentext_means %>%
  ggplot(aes(x = mean, y = variable)) + 
  geom_col(aes(fill = spending_type), width = 0.5) + 
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0.1) + 
  theme_minimal() + 
  theme(legend.position = "bottom", 
        axis.text.y = element_text(size = 11),
        legend.text = element_text(size = 11)) + 
  labs(y = NULL, x = "Share of respondents", fill = "Type of spending") +
  scale_fill_manual(values = c("#00BA38","#619CFF"), breaks = c("Relief", "Prevention"))
```

\pagebreak

## Tables

The code below generates Table 1 in the manuscript.

```{r opentext-lm-1, results = 'hide'}

# The regressions below estimate linear relationship between mentions of a concept (in open text responses) and policy preferences (from closed-ended responses).

# Regress support for spending more on relief on whether the respondent mentioned relief in open text

opentext_lm_more_relief <- lm(more_relief ~ mention_relief, data = opentext_responses)
opentext_se_more_relief <- sqrt(diag(sandwich::vcovHC(opentext_lm_more_relief, type = "HC1")))

# Regress support for spending more on prevention on whether the respondent mentioned prevention in open text

opentext_lm_more_prevention <- lm(more_prevention ~ mention_prev, data = opentext_responses)
opentext_se_more_prevention <- sqrt(diag(sandwich::vcovHC(opentext_lm_more_prevention, type = "HC1")))

# Regress support for prioritizing prevention over relief on whether the respondent mentioned relief in open text

opentext_lm_prefer_prevention_mention_relief <- lm(prefer_prevention ~ mention_relief, data = opentext_responses)
opentext_se_prefer_prevention_mention_relief <- sqrt(diag(sandwich::vcovHC(opentext_lm_prefer_prevention_mention_relief, type = "HC1")))

# Regress support for prioritizing prevention over relief on whether the respondent mentioned prevention in open text

opentext_lm_prefer_prevention_mention_prevention <- lm(prefer_prevention ~ mention_prev, data = opentext_responses)
opentext_se_prefer_prevention_mention_prevention <- sqrt(diag(sandwich::vcovHC(opentext_lm_prefer_prevention_mention_prevention, type = "HC1")))

# Combine the regression model objects and vectors containing robust standard errors separately into lists

opentext_lm_1 <- list(opentext_lm_more_relief,
                      opentext_lm_more_prevention,
                      opentext_lm_prefer_prevention_mention_relief,
                      opentext_lm_prefer_prevention_mention_prevention
                      )

opentext_se_1 <- list(opentext_se_more_relief,
                      opentext_se_more_prevention,
                      opentext_se_prefer_prevention_mention_relief,
                      opentext_se_prefer_prevention_mention_prevention
                      )

# Table 1

stargazer::stargazer(opentext_lm_1,
                     se = opentext_se_1,
                     dep.var.caption = "",
                     dep.var.labels = c("More on relief", "More on prevention", "Prioritize prevention"),
                     covariate.labels = c("Mention relief", "Mention prevention"),
                     label = "opentext-validation",
                     out = "tables/opentext-validation.tex",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     df = FALSE,
                     header = FALSE,
                     omit.stat = c("f"))
```

\input{tables/opentext-validation.tex}

## Additional Statistics

The code below shows the share of respondents who support prioritizing prevention over relief, as reported in the manuscript (see p.10).

```{r share-prefer-prevention}

survey_1 %>%
  filter(!is.na(prefer_prevention_orig)) %>%
  select(prefer_prevention_orig, weights) %>%
  mutate(yes = ifelse(prefer_prevention_orig == 1, 1, 0),
         no = ifelse(prefer_prevention_orig == 0, 1, 0),
         other = ifelse(prefer_prevention_orig == 0.5, 1, 0)
         ) %>%
  summarize(across(.cols = -c(prefer_prevention_orig, weights), 
                   .fns = ~ weighted.mean(.x, w = weights, na.rm = TRUE))) %>%
  ungroup() %>%
  mutate(across(.cols = everything(),
                .fns = ~ round(.x, digits = 2))
            ) %>%
  kable() %>%
  kable_styling(latex_options = c("scale_down", "hold_position"))
```

The code below generates the statistics on the intensity of preferences discussed in the manuscript (see p.10)

```{r overall-intensity}

# Calculate share of respondents who said they want "much more" spending on disater relief and prevention

survey_1 %>%
  select(ResponseId, weights, more_relief, more_prevention) %>%   
  mutate(across(.cols = -c(ResponseId, weights), .fns = ~ case_when(.x == 1 ~ 1, 
                                                                    is.na(.x) ~ NA_real_,
                                                                    TRUE ~ 0))) %>%
  summarize(across(.cols = -c(ResponseId, weights), .fns = ~ 100*round(weighted.mean(.x, w = weights, na.rm = TRUE), digits = 3))) %>%
  kable() %>%
  kable_styling(latex_options = c("scale_down", "hold_position"))
```

The code below generates the statistics on the intensity of preferences in the candidate evaluation discussed in the manuscript (see p.19).

```{r cand-eval-intensity}

# Recode the variables to be 1 only when the variable is equal to 1 (meaning the respondents indicated he/she would definitely vote for the candidate); then calculate averages

survey_1 %>%
  select(ResponseId, weights, cand_more_prev, cand_more_relief, cand_not_prev, cand_not_relief) %>%   
  mutate(across(.cols = -c(ResponseId, weights), .fns = ~ case_when(.x == 1 ~ 1, 
                                                                    is.na(.x) ~ NA_real_,
                                                                    TRUE ~ 0))) %>%
  summarize(across(.cols = -c(ResponseId, weights), .fns = ~ 100*round(weighted.mean(.x, w = weights, na.rm = TRUE), digits = 4))) %>%
  kable() %>%
  kable_styling(latex_options = c("scale_down", "hold_position"))
```

The code below estimates the share of respondents who mentioned natural disasters and public health disasters in their open text responses, discussed in the manuscript (p. 23) and the online appendix (Section A.7).

```{r health-natural-mentions-calc}

# Write functions for detecting each instance of a word indicating a natural or public health disaster

health_labels <- function(tbl){
  tbl <- tbl %>%
    mutate(health = str_detect(open_text, pattern = "health"),
           covid = str_detect(open_text, pattern = "covid"),
           pandemic = str_detect(open_text, pattern = "pandemic"),
           virus = str_detect(open_text, pattern = "virus"),
           vaccine = str_detect(open_text, pattern = "vaccin"),
           flu = str_detect(open_text, pattern = "flu"),
           cdc = str_detect(open_text, pattern = "cdc"),
           opioid = str_detect(open_text, pattern = "opioid"),)
  return(tbl)
  }

natural_labels <- function(tbl){
  tbl <- tbl %>%
    mutate(natural = str_detect(open_text, pattern = "natural"),
           fire = str_detect(open_text, pattern = "fire"),
           flood = str_detect(open_text, pattern = "flood"),
           earthquake = str_detect(open_text, pattern = "earthquake"),
           storm = str_detect(open_text, pattern = "storm"),
           hurricane = str_detect(open_text, pattern = "hurricane"),
           tornado = str_detect(open_text, pattern = "tornado"),
           levee = str_detect(open_text, pattern = "levee"),
           weather = str_detect(open_text, pattern = "weather"),
           seismic = str_detect(open_text, pattern = "seismic"),
           drought = str_detect(open_text, pattern = "drought"),
           warming = str_detect(open_text, pattern = "warming"),
           climate = str_detect(open_text, pattern = "climate"))
  return(tbl)
  }

# Execute the functions on the first open-ended text response

health_mentions_1 <- health_labels(open_text_data_1)
natural_mentions_1 <- natural_labels(open_text_data_1)

# Execute the functions on the second open-ended text response

health_mentions_2 <- health_labels(open_text_data_2)
natural_mentions_2 <- natural_labels(open_text_data_2)

# Create a binary variable indicating if any of the words were mentioned
health_binary_1 <- health_mentions_1 %>%
  select(-open_text) %>%
  pivot_longer(cols = -ResponseId) %>%
  group_by(ResponseId) %>%
  summarize(total = sum(value)) %>%
  ungroup() %>%
  mutate(health = ifelse(total >=1, 1, 0)) %>%
  select(ResponseId, health)

natural_binary_1 <- natural_mentions_1 %>%
  select(-open_text) %>%
  pivot_longer(cols = -ResponseId) %>%
  group_by(ResponseId) %>%
  summarize(total = sum(value)) %>%
  ungroup() %>%
  mutate(natural = ifelse(total >=1, 1, 0)) %>%
  select(ResponseId, natural)

open_text_results_1 <- health_binary_1 %>%
  left_join(natural_binary_1, by = "ResponseId") 

# Repeat the analysis for the second open-text response

health_binary_2 <- health_mentions_2 %>%
  select(-open_text) %>%
  pivot_longer(cols = -ResponseId) %>%
  group_by(ResponseId) %>%
  summarize(total = sum(value)) %>%
  ungroup() %>%
  mutate(health = ifelse(total >=1, 1, 0)) %>%
  select(ResponseId, health)

natural_binary_2 <- natural_mentions_2 %>%
  select(-open_text) %>%
  pivot_longer(cols = -ResponseId) %>%
  group_by(ResponseId) %>%
  summarize(total = sum(value)) %>%
  ungroup() %>%
  mutate(natural = ifelse(total >=1, 1, 0)) %>%
  select(ResponseId, natural)
```

The code below displays the results for the first open-ended question (before introducing prevention v relief).

```{r health-natural-mentions-1}

health_binary_1 %>%
  left_join(natural_binary_1, by = "ResponseId") %>%
  mutate(both = ifelse(health == 1 & natural == 1, 1, 0),
         neither = ifelse(health == 0 & natural == 0, 1, 0)) %>%
  summarize(across(.cols = -ResponseId, .fns = ~ round(100*mean(.x, na.rm = TRUE), digits = 0))) %>%
  kable() %>%
  kable_styling(latex_options = c("scale_down", "hold_position"))
```

The code below displays the results for the second open-ended question (after introducing prevention v relief).

```{r health-natural-mentions-2}

health_binary_2 %>%
  left_join(natural_binary_2, by = "ResponseId") %>%
  mutate(both = ifelse(health == 1 & natural == 1, 1, 0),
         neither = ifelse(health == 0 & natural == 0, 1, 0)) %>%
  summarize(across(.cols = -ResponseId, .fns = ~ round(100*mean(.x, na.rm = TRUE), digits = 0))) %>%
  kable() %>%
  kable_styling(latex_options = c("scale_down", "hold_position"))
```

The code below shows the number of times each concept (prevention and relief) is mentioned among those who mentioned it, as a pure count and as a share of the total number of words, which is mentioned in the manuscript (p. 21).

```{r mentions-thoughts}


# In hand-coded sample (no missing data for number of prevention related thoughts), among those who mention prevention, calculate the average number of thoughts and thoughts per word

opentext_responses %>%
  filter(!is.na(number_prev), !is.na(word_count), mention_prev == 1) %>%
  mutate(prev_per_word = number_prev / word_count) %>%
  summarize(avg_prev = mean(number_prev, na.rm = TRUE),
            avg_prev_per_word = mean(prev_per_word, na.rm = TRUE)) %>%
  mutate(across(.cols = everything(),
                .fns = ~ round(.x, digits = 2))) %>%
  kable() %>%
  kable_styling(latex_options = c("scale_down", "hold_position"))


# In hand-coded sample (no missing data for number of relief related thoughts), among those who mention relief, calculate the average number of thoughts and thoughts per word

opentext_responses %>%
  filter(!is.na(number_relief), !is.na(word_count), mention_relief == 1) %>%
  mutate(relief_per_word = number_relief / word_count) %>%
  summarize(avg_relief = mean(number_relief, na.rm = TRUE),
            avg_relief_per_word = mean(relief_per_word, na.rm = TRUE)) %>%
  mutate(across(.cols = everything(),
                .fns = ~ round(.x, digits = 2))) %>%
  kable() %>%
  kable_styling(latex_options = c("scale_down", "hold_position"))

```


The code below shows the share of respondents who mentioned neither prevention nor relief in their open text responses (45), along with the share that mentioned both (10), which is discussed in the manuscript (see fn. 14 on p.21 and online appendix p.15).

```{r mention-neither-both}

opentext_responses %>%
  select(ResponseId, mention_prev, mention_relief) %>%
  mutate(both = ifelse(mention_prev == 1 & mention_relief == 1, 1, 0),
         neither = ifelse(mention_prev == 0 & mention_relief == 0, 1, 0)) %>%
  summarize(across(.cols = -ResponseId, .fns = ~ round(100*mean(.x, na.rm = TRUE), digits = 0))) %>%
  kable() %>%
  kable_styling(latex_options = c("scale_down", "hold_position"))
```

The code below loads and cleans the NSCL ballot measure data (publicly available [here](https://ippsr.msu.edu/public-policy/correlates-state-policy/ncsl-ballot-measures-correlates) and included in the replication file) and defines a vector containing the IDs for the prevention-related ballot measures discussed in the manuscript (see pp. 24-25 and online appendix Section A.5, p. 19)

Note that the NCSL data ends in 2016 and does not have CA Prop 30 from 2022 (an initiative), LA Amendment 3 from 2020, or TX Prop 8 from 2019.

```{r load-ncsl-data}

# Read in the NCSL ballot measures data

ballots <- haven::read_dta("ncslballotmeasures.dta", col_select = c("ballotid", "ballotname", "st", "year", "type", "electiontype", "passed", "topicarea", "bond_meas", "budgets", "tax_rev")) %>%
  mutate(across(.cols = c(ballotid, st, year, ballotname, type, electiontype), .fns = str_trim),
         year = as.numeric(year),
         type = factor(type)) %>%
  filter(year >= 1990) 


# Define a vector containing the ballot id numbers for the disaster prevention ballot measures

disaster_ballotids <- c("CA1006", "CA1156", "LA68", "MO143", "NJ86", "OR779", "OR780", "PA5", "RI81")
```

The code below establishes that all 9 of the prevention-related ballot measures are general election referenda and are in NCSL topic areas for bond measures or budgets.

```{r characterize-prevention-measures}

ballots %>%
  filter(ballotid %in% c(disaster_ballotids)) %>%
  select(ballotid, type, electiontype, topicarea) %>%
  kable() %>%
  kable_styling(latex_options = c("scale_down", "hold_position"))
```

The code below shows that general election legislative referenda, from 1990-2016, pass at a rate of 74.7 percent.

```{r analyze-nationwide-referenda}

ballots %>%
  group_by(type, electiontype) %>%
  summarize(passage_rate = round(100*mean(passed, na.rm = TRUE), digits = 2), n = n(), ) %>%
  ungroup() %>%
  filter(type == "Legislative Referendum") %>%
  kable() %>%
  kable_styling(latex_options = c("scale_down", "hold_position"))
```

The code below constructs a comparison group for each measure by identifying bond, budget, or tax referenda from the same state, year, and election.  Note, we are only able to find comparable referenda for 7 of the 9 disaster measures.  The comparison group contains 31 measures.   

```{r identify-comparable-referenda}

# This for loop will construct a "control" group for each measure

for (j in disaster_ballotids){
  
  treatment <- ballots %>% filter(ballotid == j)

  control <- ballots %>%
    semi_join(treatment, by = c("st", "type", "year", "electiontype")) %>% # Limit to same state, ballot measure type, and election
    filter(bond_meas == 1 | budgets == 1 | tax_rev == 1, # Limit to bond, budget, or tax measures
           ballotid != j) %>% # Remove the treatment measure
    mutate(treatment_ballotid = j)

  if (j == disaster_ballotids[1]){control_referenda <- control} else {control_referenda <- control_referenda %>% bind_rows(control)}
}

# Identify the disaster referenda for which we were able to find matches and limit the "treatment" group to this set of referenda
matched_referenda_ids <- control_referenda %>% select(treatment_ballotid) %>% distinct()
matched_referenda <- ballots %>% semi_join(matched_referenda_ids, by = c("ballotid" = "treatment_ballotid"))

# Join the datsets together

combined_referenda <- control_referenda %>% 
  select(-treatment_ballotid) %>%
  mutate(prevention = 0) %>%
  bind_rows(matched_referenda %>% mutate(prevention = 1))
```

The code below illustrates the difference in mean passage rates between the prevention referenda (71.4 percent), and non-prevention spending referenda (64.5 percent).

```{r ballot-test}

lm_robust(passed ~ prevention, se_type = "stata", data = combined_referenda) %>% 
  tidy() %>%
  mutate(estimate = round(100*estimate, digits = 1)) %>%
  select(term, estimate) %>%
  kable() %>%
  kable_styling(latex_options = c("scale_down", "hold_position"))
```

\pagebreak

# Online Appendix

## Figures

The code below estimates the regression models that are plotted in Figure A.1 in the online appendix.

```{r, main-reg-lm-robust}

# For plotting, the dependent variables and demographic controls are the same.  But we do not include income or exposure in the predictors of interest (these results are reported in the raw regression output but are not discussed in connection with the figure).

main_dep_vars <- c("more_relief", "more_prevention", "prefer_prevention", "more_pandemic_relief", "more_pandemic_prevention", "prefer_pandemic_prevention")
main_demog_vars <- "age + race + gender" 
main_other_vars_plot <- c("educ", "party", "ideology", "biden", "trust_in_gov", "trust_in_media", "need_experts")

# This function will run a linear model (with robust standard errors) with demographic vars only, then add in other variables one at a time

lm_one_dep_var_robust <- function(dep_var, demog_vars, other_vars, dataset){
  models <- list()
  lm_spec <- formula(paste(dep_var, " ~ ", demog_vars, sep = ""))
  models[[1]] <- lm_robust(lm_spec, data = dataset, weights = weights, se_type = "stata")
  for (j in 1:length(other_vars)){
    lm_spec <- formula(paste(dep_var, " ~ ", demog_vars, " + ", other_vars[j], sep = ""))
    models[[j+1]] <- lm_robust(lm_spec, data = dataset, weights = weights, se_type = "stata")
  }
  return(models)
}

# This function will execute lm_one_dep_var_robust across all dependent variables (dep_vars)

lm_all_dep_vars_robust <- function(dep_vars, demog_vars, other_vars, dataset){
  all_models <- list()
  for (i in 1:length(dep_vars)){
    all_models[[i]] <- lm_one_dep_var_robust(dep_var = dep_vars[i], demog_vars = demog_vars, other_vars = other_vars, dataset = dataset)
  }
  return(all_models)
}

# Execute the functions to extract all regression output using lm_robust

main_models_robust <- lm_all_dep_vars_robust(dep_vars = main_dep_vars, demog_vars = main_demog_vars, other_vars = main_other_vars_plot, dataset = survey_1)
```

The code below extracts the estimated coefficients, standard errors, and other output (test statistics, p-values, and confidence intervals) from the `lm_robust` regressions for plotting.

```{r extract-main-reg-coefs}

# This function will extract the coefficients from a set of models for a single outcome (for plotting), while being careful to keep the estimates for demographic variables from the first model.

extract_coefs <- function(models){
  
  demog_coefs <- models[[1]] %>% tidy() # Take the first model (with demographics only) and tidy it
  demog_terms <- demog_coefs$term # Store the variable names from this model so we can exclude them from subsequent models (where we only care about the new terms)
  
  # For each remaining model, keep the estimate from the new term (taking care to also exclude party ID when the new term is exposure)
  
  for (k in 2:length(models)){ 
    
    new_coefs <- models[[k]] %>% 
      tidy() %>%
      filter(!term %in% demog_terms) %>%
      mutate(bivariate_check = as.numeric(str_detect(term, "exposure"))) # The exposure regression is conditional on party
    
    if (sum(new_coefs$bivariate_check) > 0){new_coefs <- new_coefs %>% filter(!str_detect(term, "party"))} # If the new coefs have exposure estimates, remove party estimates
    
    # Bind the new terms to the old terms
    
    if (k == 2){coefs <- demog_coefs %>% bind_rows(new_coefs)} else {coefs <- coefs %>% bind_rows(new_coefs)}
  }
  
  return(coefs)
}

# Execute the function on each set of models for each of six dependent variables

more_relief_coefs <- extract_coefs(main_models_robust[[1]])
more_prev_coefs <- extract_coefs(main_models_robust[[2]])
prioritize_prev_coefs <- extract_coefs(main_models_robust[[3]])

more_pand_relief_coefs <- extract_coefs(main_models_robust[[4]])
more_pand_prev_coefs <- extract_coefs(main_models_robust[[5]])
prioritize_pand_prev_coefs <- extract_coefs(main_models_robust[[6]])
```

The code below generates Figure A.1 in the online appendix.

```{r subgroup-analysis}

# Construct a vector of variables names for which we will exclude all estimates (at all levels)

exclude_vars <- c("(Intercept)", "race", "gender", "age")

# Bind the estimates for 3 of the 6 main outcomes together, drop the demographics, and clean the variable names (keeping the levels as is).  Use the colon (dividing variable names from levels) to separately identify the variable name for ordering (`term_group`).  Then let `term` indicate the levels.

combined_coefs <- more_prev_coefs %>%
  mutate(model = "Spend more \non prevention") %>%
  bind_rows(more_relief_coefs %>% mutate(model = "Spend more \non relief")) %>%
  bind_rows(prioritize_prev_coefs %>% mutate(model = "Prioritize \nprevention")) %>%
  group_by(model, term) %>% # Drop all levels of all excluded variables since we do not plot demographics
  filter(sum(str_detect(term, exclude_vars)) == 0) %>%
  ungroup() %>%
  group_by(term) %>%
  mutate(term = str_replace_all(term, "educ", "Education: "),
         term = str_replace_all(term, "party", "Party: "),
         term = str_replace_all(term, "ideology", "Ideology: "),
         term = str_replace_all(term, "biden", "Vote in 2020: "),
         term = str_replace_all(term, "trust_in_gov", "Trust in Gov: "),
         term = str_replace_all(term, "trust_in_media", "Trust in Sch. & Media: "),
         term = str_replace_all(term, "need_experts", "Need Experts: "),
         term_group = str_sub(term, start = 1, end = (str_locate(term, ":")[1] - 1)),
         term = str_sub(term, start = (str_locate(term, ":")[1] + 2)),
         term = case_when(term == "Bachelor's degree or higher" ~ "College deg.",
                          TRUE ~ term)) %>% 
  ungroup() %>%
  select(model, term_group, term, estimate, conf.low, conf.high)

# We want to order the figure panels so that they are in approximately increasing order by effect size according to one of the models.  Here, we determine the maximum effect size for each variable (`term_group`) in the group of models where the dependent variable is support for spending more on prevention.  We use this variable to arrange the panels.

panel_order <- combined_coefs %>%
  filter(model == "Spend more \non prevention") %>%
  group_by(term_group) %>%
  summarize(max_effect = max(estimate)) %>%
  ungroup() %>%
  arrange(max_effect) %>%
  select(term_group, max_effect) %>%
  mutate(panel_order = 1:n())

# Take the combined coefficients, convert the variable specifying the dependent variable (`model`) into a factor, and do the same for the variable names (`term_group`) using the factor ordering identified above (in `panel_order`).

ordered_coefs <- combined_coefs %>%
  mutate(model = factor(model, levels = c("Prioritize \nprevention","Spend more \non prevention","Spend more \non relief")),
         term_group = factor(term_group, levels = panel_order$term_group)) %>%
  select(model, term_group, term, estimate, conf.low, conf.high) 

# Plot the results

ordered_coefs %>%
  mutate(term_group = factor(term_group)) %>%
  ggplot(aes(x = estimate, y = term, color = model, xmin = conf.low, xmax = conf.high)) + 
  geom_point(position = position_dodge(width = 0.5)) + 
  geom_errorbar(width = 0.2, position = position_dodge(width = 0.5)) + 
  theme_minimal() + 
  geom_vline(aes(xintercept = 0), color = "blue") + 
  facet_wrap(~term_group, scales = "free_y", dir = "v") + 
  labs(x = NULL, y = NULL, color = "Outcome") +
  scale_color_manual(values = c("#00BA38","#619CFF","#F8766D"), breaks = c("Spend more \non relief", "Spend more \non prevention", "Prioritize \nprevention")) + 
  theme(legend.position = "bottom")
```

\pagebreak

The code below generates Figure A.2 in the online appendix.

```{r pandemic-subgroup-analysis}

# Bind the estimates for the remaining 3 outcomes together, drop the demographics, and clean the variable names (keeping the levels as is).  Use the colon (dividing variable names from levels) to separately identify the variable name for ordering (`term_group`).  Then let `term` indicate the levels.

combined_coefs <- more_pand_prev_coefs %>%
  mutate(model = "More on \npand. prev.") %>%
  bind_rows(more_pand_relief_coefs %>% mutate(model = "More on \npand. relief")) %>%
  bind_rows(prioritize_pand_prev_coefs %>% mutate(model = "Prioritize \npand. prev.")) %>%
  group_by(model, term) %>%
  filter(sum(str_detect(term, exclude_vars)) == 0) %>%
  ungroup() %>%
  group_by(term) %>%
  mutate(term = str_replace_all(term, "educ", "Education: "),
         term = str_replace_all(term, "party", "Party: "),
         term = str_replace_all(term, "ideology", "Ideology: "),
         term = str_replace_all(term, "biden", "Vote in 2020: "),
         term = str_replace_all(term, "trust_in_gov", "Trust in Gov: "),
         term = str_replace_all(term, "trust_in_media", "Trust in Sch. & Media: "),
         term = str_replace_all(term, "need_experts", "Need Experts: "),
         term_group = str_sub(term, start = 1, end = (str_locate(term, ":")[1] - 1)),
         term = str_sub(term, start = (str_locate(term, ":")[1] + 2)),
         term = case_when(term == "Bachelor's degree or higher" ~ "College deg.",
                          TRUE ~ term)) %>% 
  ungroup() %>%
  select(model, term_group, term, estimate, conf.low, conf.high)

# We want to order the figure panels so that they are in approximately increasing order by effect size according to one of the models.  Here, we determine the maximum effect size for each variable (`term_group`) in the group of models where the dependent variable is support for spending more on pandemic prevention.  We use this variable to arrange the panels.

panel_order <- combined_coefs %>%
  filter(model == "More on \npand. prev.") %>%
  group_by(term_group) %>%
  summarize(max_effect = max(estimate)) %>%
  ungroup() %>%
  arrange(max_effect) %>%
  select(term_group, max_effect) %>%
  mutate(panel_order = 1:n())

# Take the combined coefficients, convert the variable specifying the dependent variable (`model`) into a factor, and do the same for the variable names (`term_group`) using the factor ordering identified above (in `panel_order`).

ordered_coefs <- combined_coefs %>%
  mutate(model = factor(model, levels = c("Prioritize \npand. prev.","More on \npand. prev.","More on \npand. relief")),
         term_group = factor(term_group, levels = panel_order$term_group)) %>% 
  select(model, term_group, term, estimate, conf.low, conf.high) 

# Plot the results

ordered_coefs %>%
  mutate(term_group = factor(term_group)) %>%
  ggplot(aes(x = estimate, y = term, color = model)) + 
  geom_point(position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0.2, position = position_dodge(width = 0.5)) + 
  theme_minimal() + 
  geom_vline(aes(xintercept = 0), color = "black") + 
  facet_wrap(~term_group, scales = "free_y", dir = "v") + 
  labs(x = NULL, y = NULL, color = "Outcome") +
  scale_color_manual(values = c("#00BA38","#619CFF","#F8766D"), breaks = c("More on \npand. relief", "More on \npand. prev.", "Prioritize \npand. prev.")) + 
  theme(legend.position = "bottom")
```

\pagebreak

The code below generates Figure A.3 in the online appendix.

```{r marginal-means-pandemic}

mm_with_se %>%
  filter(dep_var %in% c("Spend more \npand. relief", "Spend more \npand. prev.","Prioritize \npand. prev.", "Spend more \nCovid prep.")) %>%
  mutate(dep_var = factor(dep_var, levels = c("Spend more \nCovid prep.","Prioritize \npand. prev.", "Spend more \npand. prev.","Spend more \npand. relief")),
         x_min = mean - 1.386*SE,
         x_max = mean + 1.386*SE) %>%
  ggplot(aes(x = mean, y = sub_group, fill = dep_var, xmin = x_min, xmax = x_max)) +
  geom_bar(width = 0.5, stat = "identity", position = position_dodge(c(0.5, 0.5, 0.9, 0.6))) + #position = "dodge", 
  geom_errorbar(width = 0.1, position = position_dodge(c(0.5, 0.5, 0.9, 0.6))) + 
  theme_minimal() +
  labs(x = NULL, y = NULL, fill = "Outcome", size = 11) +
  xlim(c(0,1)) +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  geom_text(data = annotations, aes(x = x, y = y, label = label, fill = NULL, xmin = NULL, xmax = NULL), size = 4) +
  theme(legend.position = "bottom", 
        axis.text.y = element_text(size = 11), 
        legend.text = element_text(size = 9)) + 
  scale_fill_manual(values = c("#00BA38","#619CFF","#F8766D","#C77CFF"), breaks = c("Spend more \npand. relief", "Spend more \npand. prev.","Prioritize \npand. prev.", "Spend more \nCovid prep."))
```

\pagebreak

The code below generates Figure A.4 in the online appendix.

```{r reasons-ttest}

# Create new dummy variables for whether relief or prevention went first (reversing the contrasts gives the estimated coefficients the same direction)

stability_ttest_data <- survey_1 %>% 
  select(after_reasons_prefer_prevention, after_reasons_prefer_relief, weights, reasons_prev_first) %>% 
  filter(!is.na(reasons_prev_first)) %>%
  mutate(relief_first = case_when(reasons_prev_first == 0 ~ "yes", # Create binary versions for order condition
                                  TRUE ~ "no"),
         prevention_first = case_when(reasons_prev_first == 1 ~ "yes",
                                      TRUE ~ "no"),
         across(.cols = c(prevention_first, relief_first),
                .fns = ~ factor(.x, levels = c("yes","no"))))

# Figure A.4

lm_robust(after_reasons_prefer_prevention ~ prevention_first, data = stability_ttest_data, weights = weights, se_type = "stata") %>% tidy() %>% 
  bind_rows(lm_robust(after_reasons_prefer_relief ~ relief_first, data = stability_ttest_data, weights = weights, se_type = "stata") %>% tidy()) %>%
  filter(term != "(Intercept)") %>%
  mutate(variable = case_when(outcome == "after_reasons_prefer_prevention" ~ "Prioritize prevention over relief",
                              outcome == "after_reasons_prefer_relief" ~ "Prioritize relief over prevention",
                              TRUE ~ NA_character_)) %>%
  ggplot(aes(x = estimate, y = reorder(variable, desc(variable)), xmin = conf.low, xmax = conf.high)) +
  geom_point() + 
  geom_errorbarh(height = 0.05) + 
  labs(x = NULL, y = NULL) +
  theme_bw() +
  theme(legend.position = "bottom") +
  xlim(c(-0.2, 0))
```

\pagebreak

The code below generates Figure A.5 in the online appendix.

```{r over-time-transparency}

# Calculate difference in mean support over time (between 2021 and 2023 surveys) for same wording ("natural and public health disasters")

change_over_time <- svyby(~more_spending_binary + more_relief_binary + more_prevention_binary + prefer_prevention, 
                          by = ~survey, 
                          design = combined_surveys_same_wording_design, 
                          FUN = svymean, 
                          na.rm = TRUE) %>%
  as_tibble() %>%
  pivot_longer(-survey, names_to = "variable", values_to = "estimate")

# Extract standard errors

change_over_time_stderrors <- change_over_time %>%
  filter(str_detect(variable, "se.")) %>%
  rename(std.error = estimate) %>%
  mutate(variable = str_remove(variable, "se."))

# Limit to means, add standard errors back in, and clean remaining data

change_over_time <- change_over_time %>%
  filter(!str_detect(variable, "se.")) %>%
  left_join(change_over_time_stderrors, by = c("survey", "variable")) %>%
  rename(outcome = variable) %>%
  mutate(outcome = case_when(outcome == "more_spending_binary" ~ "Spend more \non disasters",
                             outcome == "more_relief_binary" ~ "Spend more \non relief",
                             outcome == "more_prevention_binary" ~ "Spend more \non prev.",
                             outcome == "prefer_prevention" ~ "Prioritize \nprev.",
                             TRUE ~ NA_character_),
         outcome = factor(outcome, levels = c("Prioritize \nprev.", "Spend more \non prev.", "Spend more \non relief", "Spend more \non disasters")),
         year = case_when(survey == "survey_1" ~ "2021",
                          survey == "survey_2" ~ "2023",
                          TRUE ~ NA_character_),
         year = factor(year, levels = c("2023", "2021")))

# Figure A.5

change_over_time %>%
  ggplot(aes(x = estimate, y = outcome, xmin = estimate - 1.386*std.error, xmax = estimate + 1.386*std.error, fill = outcome, alpha = year)) +
  geom_col(width = 0.5, position = position_dodge(width = 0.7)) +
  theme_minimal() +
  theme(legend.position = "none") +
  geom_errorbar(width = 0.1, position = position_dodge(width = 0.7)) +
  labs(x = NULL, y = NULL, fill = NULL, alpha = NULL) + 
  scale_alpha_manual(values = c(1.0, 1.0), breaks = c("2021","2023")) + 
  scale_fill_manual(values = c("darkgrey","#00BA38","#619CFF","#F8766D"), breaks = c("Spend more \non disasters", "Spend more \non relief", "Spend more \non prev.", "Prioritize \nprev.")) + 
  xlim(c(0,1)) + 
  geom_text(aes(x = 0.1, label = year, hjust = "left"), size = 3, fontface = "bold", position = position_dodge(width = 0.7))
```

\pagebreak

The code below generates Figure A.6 in the online appendix.

```{r wording-transparency}

# Calculate difference in mean support across different question wordings in 2023 survey

condition_estimates <- svyby(~more_spending_binary + more_relief_binary + more_prevention_binary + prefer_prevention, 
                            by = ~condition,
                            design = survey_2_design,
                            FUN = svymean,
                            na.rm = TRUE) %>% 
  as_tibble() %>%
  pivot_longer(-condition, names_to = "variable", values_to = "estimate")

# Extract standard errors

condition_estimates_stderrors <- condition_estimates %>%
  filter(str_detect(variable, "se.")) %>%
  rename(std.error = estimate) %>%
  mutate(variable = str_remove(variable, "se."))

# Limit to means, add standard errors back in, and clean remaining data

condition_estimates <- condition_estimates %>%
  filter(!str_detect(variable, "se.")) %>%
  left_join(condition_estimates_stderrors, by = c("condition", "variable")) %>%
  rename(outcome = variable) %>%
  mutate(outcome = case_when(outcome == "more_spending_binary" ~ "Spend more \non disasters",
                             outcome == "more_relief_binary" ~ "Spend more \non relief",
                             outcome == "more_prevention_binary" ~ "Spend more \non prev.",
                             outcome == "prefer_prevention" ~ "Prioritize \nprev.",
                             TRUE ~ NA_character_),
         outcome = factor(outcome, levels = c("Prioritize \nprev.", "Spend more \non prev.", "Spend more \non relief", "Spend more \non disasters")),
         condition = case_when(condition == "natural_disasters" ~ "Natural disasters",
                               condition == "health_disasters" ~ "Public health\ndisasters",
                               condition == "natural_and_health_disasters" ~ "Natural and public\nhealth disasters",
                               TRUE ~ NA_character_),
         condition = factor(condition, levels = c("Natural and public\nhealth disasters", "Public health\ndisasters", "Natural disasters"))
         )

# Figure A.6

condition_estimates %>%
  filter(condition != "Natural and public\nhealth disasters") %>%
  mutate(condition = as.character(condition), 
         condition = str_replace_all(condition, pattern = "\n", " "),
         condition = factor(condition, levels = c("Public health disasters", "Natural disasters"))) %>%
  ggplot(aes(x = estimate, y = outcome, xmin = estimate - 1.386*std.error, xmax = estimate + 1.386*std.error, fill = outcome, alpha = condition)) + # 
  geom_col(width = 0.5, position = position_dodge(width = 0.7)) +
  theme_minimal() +
  theme(legend.position = "none") +
  geom_errorbar(width = 0.1, position = position_dodge(width = 0.7)) +
  labs(x = NULL, y = NULL, alpha = NULL) + # 
  scale_alpha_manual(values = c(1.0, 1.0), breaks = c("Natural disasters", "Public health disasters")) + # "Natural and public\nhealth disasters"
  scale_fill_manual(values = c("grey","#00BA38","#619CFF","#F8766D"), breaks = c("Spend more \non disasters", "Spend more \non relief", "Spend more \non prev.", "Prioritize \nprev.")) + 
  #guides(fill = "none", alpha = "none") + 
  xlim(c(0,1)) + 
  geom_text(aes(x = 0.1, label = condition, hjust = "left"), size = 3, fontface = "bold", position = position_dodge(width = 0.7))
```


## Tables

\setcounter{table}{4}
\renewcommand{\thetable}{A.\arabic{table}}

Note, Tables A.1-A.3 in the online appendix are qualitative tables that were manually assembled.  Table A.1 is a summary of prior polling data from polling databases.  Table A.2 contains variable descriptions, codings, and summary statistics for our survey.  Table A.3 is a summary of the prevention-related ballot measures we identified through searches of the NCSL database, as described in the manuscript and online appendix Section A.5.

The five code segments below reproduce Table A.4 in the online appendix. First, we gather regression results for binary outcomes with controls.

```{r binary-outcomes-with-controls}

h15 <- lm_robust(more_prevention_binary ~ condition + age + income + gender + race + educ, se_type = "stata", data = survey_2_nd_or_hd, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H15")

h16 <- lm_robust(prefer_prevention ~ condition + age + income + gender + race + educ, se_type = "stata", data = survey_2_nd_or_hd, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H16")

h17 <- lm_robust(more_prevention_binary ~ survey + age + income + gender + race + educ, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H17")

h18 <- lm_robust(prefer_prevention ~ survey + age + income + gender + race + educ, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H18")

h19 <- lm_robust(more_spending_binary ~ condition + age + income + gender + race + educ, se_type = "stata", data = survey_2, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H19")

h20 <- lm_robust(more_spending_binary ~ survey + age + income + gender + race + educ, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H20")

h21 <- lm_robust(more_relief_binary ~ condition + age + income + gender + race + educ, se_type = "stata", data = survey_2_nd_or_hd, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H21")

h22 <- lm_robust(more_relief_binary ~ survey + age + income + gender + race + educ, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H22")

binary_with_controls <- h15 %>%
  bind_rows(h16, h17, h18, h19, h20, h21, h22) %>%
  mutate(across(.cols = c(estimate, std.error), 
                .fns = ~ round(.x, digits = 3)),
         estimate = case_when((p.value > 0.01 & p.value <= 0.05) ~ paste(estimate, "*", "\n(", std.error, ")", sep = ""),
                              (p.value <= 0.01) ~ paste(estimate,  "**", "\n(", std.error, ")", sep = ""),
                              (p.value > 0.05) ~ paste(estimate,  "\n(", std.error, ")", sep = ""),
                              TRUE ~ NA_character_)
         ) %>%
  rename(`Binary, controls` = estimate) %>%
  select(term, outcome, Hypothesis, `Binary, controls`)
```

Next, we gather regression results for binary outcomes with no controls.

```{r binary-outcomes-with-no-controls}

h15 <- lm_robust(more_prevention_binary ~ condition, se_type = "stata", data = survey_2_nd_or_hd, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H15")

h16 <- lm_robust(prefer_prevention ~ condition, se_type = "stata", data = survey_2_nd_or_hd, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H16")

h17 <- lm_robust(more_prevention_binary ~ survey, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H17")

h18 <- lm_robust(prefer_prevention ~ survey, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H18")

h19 <- lm_robust(more_spending_binary ~ condition, se_type = "stata", data = survey_2, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H19")

h20 <- lm_robust(more_spending_binary ~ survey, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H20")

h21 <- lm_robust(more_relief_binary ~ condition, se_type = "stata", data = survey_2_nd_or_hd, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H21")

h22 <- lm_robust(more_relief_binary ~ survey, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H22")

binary_with_no_controls <- h15 %>%
  bind_rows(h16, h17, h18, h19, h20, h21, h22) %>%
  mutate(across(.cols = c(estimate, std.error), 
                .fns = ~ round(.x, digits = 3)),
         estimate = case_when((p.value > 0.01 & p.value <= 0.05) ~ paste(estimate, "*", "\n(", std.error, ")", sep = ""),
                              (p.value <= 0.01) ~ paste(estimate,  "**", "\n(", std.error, ")", sep = ""),
                              (p.value > 0.05) ~ paste(estimate,  "\n(", std.error, ")", sep = ""),
                              TRUE ~ NA_character_)
         ) %>%
  rename(`Binary, no controls` = estimate) %>%
  select(Hypothesis, `Binary, no controls`)
```

Then, we gather regression results for outcomes on the five-point scale with controls.

```{r five-point-with-controls}

# Note, H16 and H18, for prioritizing prevention over relief, is only binary and does not have a five-point measure.

h15 <- lm_robust(more_prevention ~ condition + age + income + gender + race + educ, se_type = "stata", data = survey_2_nd_or_hd, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H15")

h17 <- lm_robust(more_prevention ~ survey + age + income + gender + race + educ, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H17")

h19 <- lm_robust(more_spending ~ condition + age + income + gender + race + educ, se_type = "stata", data = survey_2, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H19")

h20 <- lm_robust(more_spending ~ survey + age + income + gender + race + educ, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H20")

h21 <- lm_robust(more_relief ~ condition + age + income + gender + race + educ, se_type = "stata", data = survey_2_nd_or_hd, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H21")

h22 <- lm_robust(more_relief ~ survey + age + income + gender + race + educ, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H22")

five_point_with_controls <- h15 %>%
  bind_rows(h17, h19, h20, h21, h22) %>%
  mutate(across(.cols = c(estimate, std.error), 
                .fns = ~ round(.x, digits = 3)),
         estimate = case_when((p.value > 0.01 & p.value <= 0.05) ~ paste(estimate, "*", "\n(", std.error, ")", sep = ""),
                              (p.value <= 0.01) ~ paste(estimate,  "**", "\n(", std.error, ")", sep = ""),
                              (p.value > 0.05) ~ paste(estimate,  "\n(", std.error, ")", sep = ""),
                              TRUE ~ NA_character_)
         ) %>%
  rename(`Five-point, controls` = estimate) %>%
  select(Hypothesis, `Five-point, controls`)
```

Then, we gather regression results for outcomes on the five-point scale with no controls.

```{r five-point-with-no-controls}

h15 <- lm_robust(more_prevention ~ condition, se_type = "stata", data = survey_2_nd_or_hd, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H15")

h17 <- lm_robust(more_prevention ~ survey, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H17")

h19 <- lm_robust(more_spending ~ condition, se_type = "stata", data = survey_2, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H19")

h20 <- lm_robust(more_spending ~ survey, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H20")

h21 <- lm_robust(more_relief_binary ~ condition, se_type = "stata", data = survey_2_nd_or_hd, weights = weights) %>% 
  tidy() %>%
  filter(term == "conditionhealth_disasters") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H21")

h22 <- lm_robust(more_relief_binary ~ survey, se_type = "stata", data = combined_survey_same_wording, weights = weights) %>% 
  tidy() %>%
  filter(term == "surveysurvey_2") %>%
  select(term, estimate, std.error, p.value, df, outcome) %>%
  mutate(Hypothesis = "H22")

five_point_with_no_controls <- h15 %>%
  bind_rows(h17, h19, h20, h21, h22) %>%
  mutate(across(.cols = c(estimate, std.error), 
                .fns = ~ round(.x, digits = 3)),
         estimate = case_when((p.value > 0.01 & p.value <= 0.05) ~ paste(estimate, "*", "\n(", std.error, ")", sep = ""),
                              (p.value <= 0.01) ~ paste(estimate,  "**", "\n(", std.error, ")", sep = ""),
                              (p.value > 0.05) ~ paste(estimate,  "\n(", std.error, ")", sep = ""),
                              TRUE ~ NA_character_)
         ) %>%
  rename(`Five-point, no controls` = estimate) %>%
  select(Hypothesis, `Five-point, no controls`)
```

Last, we generate the summary table.

```{r publish-summary-table, warning = FALSE, message = FALSE}

binary_with_controls %>%
  left_join(binary_with_no_controls, by = c("Hypothesis")) %>%
  left_join(five_point_with_controls, by = c("Hypothesis")) %>%
  left_join(five_point_with_no_controls, by = c("Hypothesis")) %>%
  mutate(term = case_when(term == "conditionhealth_disasters" ~ "Health disasters compared to natural disasters",
                          term == "surveysurvey_2" ~ "2023 compared to 2021"),
         Outcome = case_when(outcome == "more_prevention_binary" ~ "Spend more on prev.",
                                 outcome == "prefer_prevention" ~ "Prioritize prev.",
                                 outcome == "more_relief_binary" ~ "Spend more on relief",
                                 outcome == "more_spending_binary" ~ "Spend more on disasters"
                                 ),
         term = factor(term, levels = c("2023 compared to 2021", "Health disasters compared to natural disasters")),
         Outcome = factor(Outcome, levels = c("Spend more on disasters", "Spend more on relief", "Spend more on prev.", "Prioritize prev."))
         ) %>%
  arrange(term, Outcome) %>% 
  rename(`Hyp.` = Hypothesis) %>%
  select(Outcome, `Hyp.`, `Binary, controls`, `Binary, no controls`,`Five-point, controls`, `Five-point, no controls`) %>%
  mutate(`Hyp. supported?` = "Yes",
         across(.cols = c(`Binary, controls`, `Binary, no controls`,`Five-point, controls`, `Five-point, no controls`),
                .fns = ~ linebreak(.x, align = "c")),
         across(.cols = c(`Five-point, controls`, `Five-point, no controls`),
                .fns = ~ replace_na(.x, "--"))
         )%>%
  kable(escape = FALSE, align = c("l", "c", "c", "c", "c", "c", "c")) %>%
  column_spec(1, width = "1.5in") %>%
  column_spec(2, width = "0.5in") %>%
  column_spec(3:6, width = "0.7in") %>%
  column_spec(7, width = "0.5in") %>%
  pack_rows("Support in 2023 not substantially less than in 2021", 1, 4) %>%
  pack_rows("Support for health disasters not substantially greater than natural disasters", 5, 8) %>%
  footnote(general = "* p < 0.05; ** p< 0.01", threeparttable = TRUE)
```

The code below creates new variables that specify the main dependent variables or outcomes of interest for the regression analysis presented in Tables A.5-A.10 (6 variables total), the main demographic variables to include in all regressions (3 total), and the remaining independent variables which will enter the model one at a time (9 total).  Note that, in the case of exposure, we also control for partisan identification per the pre-analysis plan.

```{r set-main-reg-vars}

main_dep_vars <- c("more_relief", "more_prevention", "prefer_prevention", "more_pandemic_relief", "more_pandemic_prevention", "prefer_pandemic_prevention")
main_demog_vars <- "age + race + gender" 
main_other_vars <- c("income", "educ", "party", "ideology", "biden", "trust_in_gov", "trust_in_media", "need_experts", "party + exposure")
```

The code below estimates the regression models that are used to generate Tables A.5-A.10.  

```{r main-reg-lm}

# This function will run a linear model with demographic vars only, then add in other variables one at a time; it saves robust standard errors in the second element of a list

lm_one_dep_var <- function(dep_var, demog_vars, other_vars, dataset){
  
  models <- list()
  std_errors <- list()
  
  lm_spec <- formula(paste(dep_var, " ~ ", demog_vars, sep = ""))
  
  models[[1]] <- lm(lm_spec, data = dataset, weights = weights)
  std_errors[[1]] <- sqrt(diag(sandwich::vcovHC(models[[1]], type = "HC1")))

  for (j in 1:length(other_vars)){
    lm_spec <- formula(paste(dep_var, " ~ ", demog_vars, " + ", other_vars[j], sep = ""))
    models[[j+1]] <- lm(lm_spec, data = dataset, weights = weights)
    std_errors[[j+1]] <- sqrt(diag(sandwich::vcovHC(models[[j+1]], type = "HC1")))
    
  }
  
  output <- list()
  output[[1]] <- models
  output[[2]] <- std_errors
  
  return(output)
}

# This function will execute lm_one_dep_var across all dependent variables (dep_vars)

lm_all_dep_vars <- function(dep_vars, demog_vars, other_vars, dataset){
  all_models <- list()
  for (i in 1:length(dep_vars)){
    all_models[[i]] <- lm_one_dep_var(dep_var = dep_vars[i], demog_vars = demog_vars, other_vars = other_vars, dataset = dataset)
  }
  return(all_models)
}

# Execute the functions to extract all regression output using lm

main_models <- lm_all_dep_vars(dep_vars = main_dep_vars, demog_vars = main_demog_vars, other_vars = main_other_vars, dataset = survey_1)
```

\pagebreak

The code below generates Table A.5 in the online appendix.

```{r main-reg-tables-1, results = 'hide'}

# Create a vector of clean variable names

covariate_labels <- c("Age: 40-59", "Age: 60 or older", 
                      "Race: Black", "Race: Asian", "Race: Multiple and Other",
                      "Gender: Female", "Gender: Other", 
                      "Income: 75-149k", "Income: 150k or more",
                      "Educ: Some college", "Educ: Bachelor's or higher", 
                      "Party: Independent", "Party: Democrat", 
                      "Ideology: Moderate", "Ideology: Liberal",
                      "Voted for Other", "Did Not Vote", "Voted for Biden", 
                      "Trust in gov: Moderate", "Trust in gov: High", 
                      "Trust in media: Moderate", "Trust in media: High", 
                      "Need experts: Moderate", "Need experts: High", 
                      "Exposure: Covid", "Exposure: Disaster", "Exposure: Both")

# Table A.5

stargazer::stargazer(main_models[[1]][[1]],
                     se = main_models[[1]][[2]],
                     covariate.labels = covariate_labels,
                     dep.var.caption = "",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = c(""),
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/main-model-relief.tex",
                     label = "main-model-relief")
```

\input{tables/main-model-relief.tex}

\pagebreak

The code below generates Table A.6 in the online appendix.

```{r main-reg-tables-2, results = 'hide'}

# Table A.6

stargazer::stargazer(main_models[[2]][[1]],
                     se = main_models[[2]][[2]],
                     covariate.labels = covariate_labels,
                     dep.var.caption = "",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = c(""),
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/main-model-prevention.tex",
                     label = "main-model-prevention")
```

\input{tables/main-model-prevention.tex}

\pagebreak

The code below generates Table A.7 in the online appendix.

```{r main-reg-tables-3, results = 'hide'}

# Table A.7

stargazer::stargazer(main_models[[3]][[1]], 
                     se = main_models[[3]][[2]],
                     covariate.labels = covariate_labels,
                     dep.var.caption = "",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = c(""),
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/main-model-prioritize-prevention.tex",
                     label = "main-model-prioritize-prevention")
```

\input{tables/main-model-prioritize-prevention.tex}

\pagebreak

The code below generates Table A.8 in the appendix.

```{r main-reg-tables-4, results = 'hide'}

# Table A.8

stargazer::stargazer(main_models[[4]][[1]], 
                     se = main_models[[4]][[2]],
                     covariate.labels = covariate_labels,
                     dep.var.caption = "",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = c(""),
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/main-model-pandemic-relief.tex",
                     label = "main-model-pandemic-relief")
```

\input{tables/main-model-pandemic-relief.tex}

\pagebreak

The code below generates Table A.9 in the online appendix.

```{r main-reg-tables-5, results = 'hide'}

# Table A.9

stargazer::stargazer(main_models[[5]][[1]], 
                     se = main_models[[5]][[2]],
                     covariate.labels = covariate_labels,
                     dep.var.caption = "",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = c(""),
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/main-model-pandemic-prevention.tex",
                     label = "main-model-pandemic-prevention")
```

\input{tables/main-model-pandemic-prevention.tex}

\pagebreak

The code below generates Table A.10 in the online appendix.

```{r main-reg-tables-6, results = 'hide'}

# Table A.10

stargazer::stargazer(main_models[[6]][[1]], 
                     se = main_models[[6]][[2]],
                     covariate.labels = covariate_labels,
                     dep.var.caption = "",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = c(""),
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/main-model-prioritize-pandemic-prevention.tex",
                     label = "main-model-prioritize-pandemic-prevention")
```

\input{tables/main-model-prioritize-pandemic-prevention.tex}

\pagebreak

The code below generates the regression output that is presented in Tables A.11-A.14 of the online appendix.

```{r follow-up-regs}

# Run regressions for binary outcomes (with and without controls).

h15_binary_controls <- lm(more_prevention_binary ~ condition + age + income + gender + race + educ, data = survey_2_nd_or_hd, weights = weights) 
h15_binary_controls_ses <- sqrt(diag(sandwich::vcovHC(h15_binary_controls, type = "HC1")))

h15_binary <- lm(more_prevention_binary ~ condition, data = survey_2_nd_or_hd, weights = weights) 
h15_binary_ses <- sqrt(diag(sandwich::vcovHC(h15_binary, type = "HC1")))

h16_binary_controls <- lm(prefer_prevention ~ condition + age + income + gender + race + educ, data = survey_2_nd_or_hd, weights = weights) 
h16_binary_controls_ses <- sqrt(diag(sandwich::vcovHC(h16_binary_controls, type = "HC1")))

h16_binary <- lm(prefer_prevention ~ condition, data = survey_2_nd_or_hd, weights = weights) 
h16_binary_ses <- sqrt(diag(sandwich::vcovHC(h16_binary, type = "HC1")))

h17_binary_controls <- lm(more_prevention_binary ~ survey + age + income + gender + race + educ, data = combined_survey_same_wording, weights = weights) 
h17_binary_controls_ses <- sqrt(diag(sandwich::vcovHC(h17_binary_controls, type = "HC1")))

h17_binary <- lm(more_prevention_binary ~ survey, data = combined_survey_same_wording, weights = weights) 
h17_binary_ses <- sqrt(diag(sandwich::vcovHC(h17_binary, type = "HC1")))

h18_binary_controls <- lm(prefer_prevention ~ survey + age + income + gender + race + educ, data = combined_survey_same_wording, weights = weights) 
h18_binary_controls_ses <- sqrt(diag(sandwich::vcovHC(h18_binary_controls, type = "HC1")))

h18_binary <- lm(prefer_prevention ~ survey, data = combined_survey_same_wording, weights = weights) 
h18_binary_ses <- sqrt(diag(sandwich::vcovHC(h18_binary, type = "HC1")))

h19_binary_controls <- lm(more_spending_binary ~ condition + age + income + gender + race + educ, data = survey_2, weights = weights) 
h19_binary_controls_ses <- sqrt(diag(sandwich::vcovHC(h19_binary_controls, type = "HC1")))

h19_binary <- lm(more_spending_binary ~ condition, data = survey_2, weights = weights) 
h19_binary_ses <- sqrt(diag(sandwich::vcovHC(h19_binary, type = "HC1")))

h20_binary_controls <- lm(more_spending_binary ~ survey + age + income + gender + race + educ, data = combined_survey_same_wording, weights = weights) 
h20_binary_controls_ses <- sqrt(diag(sandwich::vcovHC(h20_binary_controls, type = "HC1")))

h20_binary <- lm(more_spending_binary ~ survey, data = combined_survey_same_wording, weights = weights) 
h20_binary_ses <- sqrt(diag(sandwich::vcovHC(h20_binary, type = "HC1")))

h21_binary_controls <- lm(more_relief_binary ~ condition + age + income + gender + race + educ, data = survey_2_nd_or_hd, weights = weights) 
h21_binary_controls_ses <- sqrt(diag(sandwich::vcovHC(h21_binary_controls, type = "HC1")))

h21_binary <- lm(more_relief_binary ~ condition, data = survey_2_nd_or_hd, weights = weights) 
h21_binary_ses <- sqrt(diag(sandwich::vcovHC(h21_binary, type = "HC1")))

h22_binary_controls <- lm(more_relief_binary ~ survey + age + income + gender + race + educ, data = combined_survey_same_wording, weights = weights)
h22_binary_controls_ses <- sqrt(diag(sandwich::vcovHC(h22_binary_controls, type = "HC1")))

h22_binary <- lm(more_relief_binary ~ survey, data = combined_survey_same_wording, weights = weights)
h22_binary_ses <- sqrt(diag(sandwich::vcovHC(h22_binary, type = "HC1")))

# Run regressions for five-point outcomes (with and without controls).

h15_five_point_controls <- lm(more_prevention ~ condition + age + income + gender + race + educ, data = survey_2_nd_or_hd, weights = weights) 
h15_five_point_controls_ses <- sqrt(diag(sandwich::vcovHC(h15_five_point_controls, type = "HC1")))

h15_five_point <- lm(more_prevention ~ condition, data = survey_2_nd_or_hd, weights = weights) 
h15_five_point_ses <- sqrt(diag(sandwich::vcovHC(h15_five_point, type = "HC1")))

h17_five_point_controls <- lm(more_prevention ~ survey + age + income + gender + race + educ, data = combined_survey_same_wording, weights = weights) 
h17_five_point_controls_ses <- sqrt(diag(sandwich::vcovHC(h17_five_point_controls, type = "HC1")))

h17_five_point <- lm(more_prevention ~ survey, data = combined_survey_same_wording, weights = weights) 
h17_five_point_ses <- sqrt(diag(sandwich::vcovHC(h17_five_point, type = "HC1")))

h19_five_point_controls <- lm(more_spending ~ condition + age + income + gender + race + educ, data = survey_2, weights = weights) 
h19_five_point_controls_ses <- sqrt(diag(sandwich::vcovHC(h19_five_point_controls, type = "HC1")))

h19_five_point <- lm(more_spending ~ condition, data = survey_2, weights = weights) 
h19_five_point_ses <- sqrt(diag(sandwich::vcovHC(h19_five_point, type = "HC1")))

h20_five_point_controls <- lm(more_spending ~ survey + age + income + gender + race + educ, data = combined_survey_same_wording, weights = weights) 
h20_five_point_controls_ses <- sqrt(diag(sandwich::vcovHC(h20_five_point_controls, type = "HC1")))

h20_five_point <- lm(more_spending ~ survey, data = combined_survey_same_wording, weights = weights) 
h20_five_point_ses <- sqrt(diag(sandwich::vcovHC(h20_five_point, type = "HC1")))

h21_five_point_controls <- lm(more_relief ~ condition + age + income + gender + race + educ, data = survey_2_nd_or_hd, weights = weights) 
h21_five_point_controls_ses <- sqrt(diag(sandwich::vcovHC(h21_five_point_controls, type = "HC1")))

h21_five_point <- lm(more_relief ~ condition, data = survey_2_nd_or_hd, weights = weights) 
h21_five_point_ses <- sqrt(diag(sandwich::vcovHC(h21_five_point, type = "HC1")))

h22_five_point_controls <- lm(more_relief ~ survey + age + income + gender + race + educ, data = combined_survey_same_wording, weights = weights)
h22_five_point_controls_ses <- sqrt(diag(sandwich::vcovHC(h22_five_point_controls, type = "HC1")))

h22_five_point <- lm(more_relief ~ survey, data = combined_survey_same_wording, weights = weights)
h22_five_point_ses <- sqrt(diag(sandwich::vcovHC(h22_five_point, type = "HC1")))
```

The code below generates Table A.11 in the online appendix.

```{r follow-up-tabs-1, results = 'hide'}

h15_and_h16 <- list(h15_binary, h15_binary_controls, h15_five_point, h15_five_point_controls,
                    h16_binary, h16_binary_controls)

h15_and_h16_ses <- list(h15_binary_ses, h15_binary_controls_ses, h15_five_point_ses, h15_five_point_controls_ses,
                        h16_binary_ses, h16_binary_controls_ses)

# Table A.11

covariate_labels <- c("Condition: Health Disasters",
                      "Age: 40-59", "Age: 60 or older",
                      "Income: 75-149k", "Income: 150k or more",
                      "Gender: Female", "Gender: Other",
                      "Race: Black", "Race: Asian", "Race: Multiple and Other",
                      "Educ: Some college", "Educ: Bachelor's or higher")

dep_var_labels <- c("More on prev. (binary)", "More on prev. (five-point scale, 0-1)", "Prioritize prev. (binary)")

stargazer::stargazer(h15_and_h16,
                     se = h15_and_h16_ses,
                     covariate.labels = covariate_labels,
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = dep_var_labels,
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/h15-and-h16.tex",
                     label = "h15-and-h16")
```

\input{tables/h15-and-h16.tex}

\pagebreak

The code below generates Table A.12 in the online appendix.

```{r follow-up-tabs-2, results = 'hide'}

h17_and_h18 <- list(h17_binary, h17_binary_controls, h17_five_point, h17_five_point_controls,
                    h18_binary, h18_binary_controls)

h17_and_h18_ses <- list(h17_binary_ses, h17_binary_controls_ses, h17_five_point_ses, h17_five_point_controls_ses,
                        h18_binary_ses, h18_binary_controls_ses)

# Table A.12

covariate_labels <- c("Survey year: 2023",
                      "Age: 40-59", "Age: 60 or older",
                      "Income: 75-149k", "Income: 150k or more",
                      "Gender: Female", "Gender: Other",
                      "Race: Black", "Race: Asian", "Race: Multiple and Other",
                      "Educ: Some college", "Educ: Bachelor's or higher")

dep_var_labels <- c("More on prev. (binary)", "More on prev. (five-point scale, 0-1)", "Prioritize prev. (binary)")

stargazer::stargazer(h17_and_h18,
                     se = h17_and_h18_ses,
                     covariate.labels = covariate_labels,
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = dep_var_labels,
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/h17-and-h18.tex",
                     label = "h17-and-h18")
```

\input{tables/h17-and-h18.tex}

\pagebreak

The code below generates Table A.13 in the online appendix.

```{r follow-up-tabs-3, results = 'hide'}

h19_and_h21 <- list(h19_binary, h19_binary_controls, h19_five_point, h19_five_point_controls,
                    h21_binary, h21_binary_controls, h21_five_point, h21_five_point_controls)

h19_and_h21_ses <- list(h19_binary_ses, h19_binary_controls_ses, h19_five_point_ses, h19_five_point_controls_ses,
                        h21_binary_ses, h21_binary_controls_ses, h21_five_point_ses, h21_five_point_controls_ses)

# Table A.13

covariate_labels <- c("Condition: Health Disasters",
                      "Age: 40-59", "Age: 60 or older",
                      "Income: 75-149k", "Income: 150k or more",
                      "Gender: Female", "Gender: Other",
                      "Race: Black", "Race: Asian", "Race: Multiple and Other",
                      "Educ: Some college", "Educ: Bachelor's or higher")

dep_var_labels <- c("More on disasters (binary)", "More on disasters (five-point scale, 0-1)", "More on relief (binary)", "More on relief (five-point scale, 0-1)")

stargazer::stargazer(h19_and_h21,
                     se = h19_and_h21_ses,
                     covariate.labels = covariate_labels,
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = dep_var_labels,
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/h19-and-h21.tex",
                     label = "h19-and-h21")
```

\input{tables/h19-and-h21.tex}

\pagebreak

The code below generate Tables A.14 in the online appendix.

```{r follow-up-tabs-4, results = 'hide'}

h20_and_h22 <- list(h20_binary, h20_binary_controls, h20_five_point, h20_five_point_controls,
                    h22_binary, h22_binary_controls, h22_five_point, h22_five_point_controls)

h20_and_h22_ses <- list(h20_binary_ses, h20_binary_controls_ses, h20_five_point_ses, h20_five_point_controls_ses,
                        h22_binary_ses, h22_binary_controls_ses, h22_five_point_ses, h22_five_point_controls_ses)

# Table A.14

covariate_labels <- c("Survey year: 2023",
                      "Age: 40-59", "Age: 60 or older",
                      "Income: 75-149k", "Income: 150k or more",
                      "Gender: Female", "Gender: Other",
                      "Race: Black", "Race: Asian", "Race: Multiple and Other",
                      "Educ: Some college", "Educ: Bachelor's or higher")

dep_var_labels <- c("More on disasters (binary)", "More on disasters (five-point scale, 0-1)", "More on relief (binary)", "More on relief (five-point scale, 0-1)")

stargazer::stargazer(h20_and_h22,
                     se = h20_and_h22_ses,
                     covariate.labels = covariate_labels,
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = dep_var_labels,
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/h20-and-h22.tex",
                     label = "h20-and-h22")
```

\input{tables/h20-and-h22.tex}

\pagebreak

The code below generates Table A.15 in the online appendix.

```{r across-reasons-types, results = 'hide'}

# Take all of the persuasiveness ratings (questions Q8_1 through Q8_9 for prevention and Q11_1 through Q11_10 for relief), pivot to long format, and create a dummy variable (spending_type) indicating whether the spending type is for prevention or relief.  Preserve the ResponseId (for clustering standard errors) and the survey weights.  

across_reason_types <- survey_1 %>%
  select(ResponseId, weights, Q8_1, Q8_2, Q8_3, Q8_4, Q8_5, Q8_6, Q8_7, Q8_8, Q8_9, Q11_1, Q11_2, Q11_3, Q11_4, Q11_5, Q11_6, Q11_7, Q11_8, Q11_9, Q11_10) %>%
  pivot_longer(-c(ResponseId, weights)) %>%
  mutate(spending_type = case_when(str_sub(name, start = 1, end = 2) == "Q8" ~ "prevention",
                           str_sub(name, start = 1, end = 3) == "Q11" ~ "relief",
                           TRUE ~ NA_character_
                           ),
         spending_type = factor(spending_type, levels = c("prevention", "relief"))) %>%
  rename(question = name, persuasiveness = value)

# There 2,828 unique ResponseIds with some data on persuasiveness (used to specify degrees of freedom for clustering below)

# across_reason_types %>% filter(!is.na(persuasiveness)) %>% select(ResponseId) %>% n_distinct()

# Regress persuasiveness on the dummy variable for spending type.

across_reason_types_lm <- lm(persuasiveness ~ spending_type, data = across_reason_types, weights = weights) 
across_reason_types_se <- lmtest::coeftest(across_reason_types_lm, 
                                           vcov = sandwich::vcovCL, 
                                           type = "HC1",
                                           df = 2827,
                                           cluster = ~ResponseId)[,2]


# Table A.15

stargazer::stargazer(across_reason_types_lm,
                     se = list(across_reason_types_se),
                     covariate.labels = c("Reason type: Relief"),
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = c("Persuasiveness of reasons"),
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/across-reason-types.tex",
                     label = "across-reason-types")
```


\input{tables/across-reason-types.tex}

<!-- \pagebreak -->

The code below generates Table A.16 in the online appendix.

```{r stability-ttest-reg, results = 'hide'}

# Regress preferences on dummy variable indicating reasons for that policy came second (baseline condition being the reasons came first)

reasons_shift_prev_lm <- lm(after_reasons_prefer_prevention ~ prevention_first, data = stability_ttest_data, weights = weights) 
reasons_shift_prev_se <- sqrt(diag(sandwich::vcovHC(reasons_shift_prev_lm, type = "HC1")))

reasons_shift_relief_lm <- lm(after_reasons_prefer_relief ~ relief_first, data = stability_ttest_data, weights = weights) 
reasons_shift_relief_se <- sqrt(diag(sandwich::vcovHC(reasons_shift_relief_lm, type = "HC1")))

stability_ttest_lm <- list(reasons_shift_prev_lm, reasons_shift_relief_lm)
stability_ttest_se <- list(reasons_shift_prev_se, reasons_shift_relief_se)

stargazer::stargazer(stability_ttest_lm, 
                     se = stability_ttest_se,
                     covariate.labels = c("Order: Prev. Second", "Order: Relief Second"),
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = c("Prioritize Prev.", "Prioritize Relief"),
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/stability-ttest.tex",
                     label = "stability-ttest")
```

\input{tables/stability-ttest.tex}

<!-- \pagebreak -->

The code below generates Table A.17 in the online appendix.

```{r reasons-reg-1, results = 'hide'}

# There are 2,459 respondents in this data set with non-missing data on candidate evaluations

# candidate_data %>% select(ResponseId) %>% n_distinct()

# Run regressions

candidate_relief_lm <- lm(support ~ cand_voted + spendingposition * spendingtype, data = candidate_relief_contrast, weights = weights)
candidate_relief_se <- lmtest::coeftest(candidate_relief_lm, vcov = sandwich::vcovCL, df = 2458, cluster = ~ResponseId, type = "HC1")[,2]

candidate_prev_lm <- lm(support ~ cand_voted + spendingposition * spendingtype, data = candidate_prev_contrast, weights = weights)
candidate_prev_se <- lmtest::coeftest(candidate_prev_lm, vcov = sandwich::vcovCL, df = 2458, cluster = ~ResponseId, type = "HC1")[,2]

# Note, the regression here is slightly different from above, because the "first" candidates only included the "spend more" on prevention / relief conditions (not the "did not spend more" conditions)
# As a result, so we cannot control for `spendingposition` (and we don't need to interact it with `spendingtype`).

candidate_first_prev_lm <- lm(support ~ cand_voted + spendingtype, data = candidate_first_prev_contrast, weights = weights)
candidate_first_prev_se <- lmtest::coeftest(candidate_first_prev_lm, vcov = sandwich::vcovHC, type = "HC1")[,2]

candidate_first_relief_lm <- lm(support ~ cand_voted + spendingtype, data = candidate_first_relief_contrast, weights = weights)
candidate_first_relief_se <- lmtest::coeftest(candidate_first_relief_lm, vcov = sandwich::vcovHC, type = "HC1")[,2]

# Print regression output

# Table A.17

stargazer::stargazer(list(candidate_relief_lm, candidate_prev_lm),
                     se = list(candidate_relief_se, candidate_prev_se),
                     covariate.labels = c("Voted for in last election: Yes", "Increased spending: Yes", 
                                          "Spending type: Prevention" , "Increased spending: Yes X Spending Type: Prevention",
                                          "Spending type: Relief" , "Increased spending: Yes X Spending Type: Relief"),
                     dep.var.caption = "",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = c("Likelihood of voting for candidate"),
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/candidate-selection.tex",
                     label = "candidate-selection.tex")
```

\input{tables/candidate-selection.tex}

\pagebreak

The code below generates Table A.18 in the online appendix.

```{r reasons-reg-2, results = 'hide'}

# Table A.18

stargazer::stargazer(candidate_first_relief_lm,
                     se = list(candidate_first_relief_se),
                     covariate.labels = c("Voted for in last election: Yes", "Type of spending increase: Prevention", 
                                          "Type of spending increase: Relief"),
                     dep.var.caption = "",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     dep.var.labels = c("Likelihood of voting for candidate"),
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/candidate-first-selection.tex",
                     label = "candidate-first-selection")

```

\input{tables/candidate-first-selection.tex}

\pagebreak

The code below generates Table A.19 in the online appendix.

```{r exposure-tests, results = 'hide'}

# Regress support for more relief, more prevention, and prioritizing prevention on demographics, party ID, and covid exposure

exposure_more_relief_covid_lm <- lm(more_relief ~ age + race + gender + party + covid_exposure, 
                                    data = survey_1,
                                    weights = weights)

exposure_more_relief_covid_se <- lmtest::coeftest(exposure_more_relief_covid_lm, 
                                                       vcov = sandwich::vcovHC, 
                                                       type = "HC1")[,2]

exposure_more_prevention_covid_lm <- lm(more_prevention ~ age + race + gender + party + covid_exposure, 
                                        data = survey_1, 
                                        weights = weights)

exposure_more_prevention_covid_se <- lmtest::coeftest(exposure_more_prevention_covid_lm, 
                                                           vcov = sandwich::vcovHC, 
                                                           type = "HC1")[,2]

exposure_prefer_prevention_covid_lm <- lm(prefer_prevention ~ age + race + gender + party + covid_exposure, 
                                          data = survey_1, 
                                          weights = weights)

exposure_prefer_prevention_covid_se <- lmtest::coeftest(exposure_prefer_prevention_covid_lm, 
                                                             vcov = sandwich::vcovHC, 
                                                             type = "HC1")[,2]

# Regress support for more relief, more prevention, and prioritizing prevention on demographics, party ID, and disaster exposure

exposure_more_relief_disaster_lm <- lm(more_relief ~ age + race + gender + party + disaster_exposure, 
                                       data = survey_1,
                                       weights = weights)

exposure_more_relief_disaster_se <- lmtest::coeftest(exposure_more_relief_disaster_lm, 
                                                          vcov = sandwich::vcovHC, 
                                                          type = "HC1")[,2]

exposure_more_prevention_disaster_lm <- lm(more_prevention ~ age + race + gender + party + disaster_exposure, 
                                           data = survey_1, 
                                           weights = weights)

exposure_more_prevention_disaster_se <- lmtest::coeftest(exposure_more_prevention_disaster_lm, 
                                                     vcov = sandwich::vcovHC, 
                                                     type = "HC1")[,2]

exposure_prefer_prevention_disaster_lm <- lm(prefer_prevention ~ age + race + gender + party + disaster_exposure, 
                                          data = survey_1, 
                                          weights = weights)

exposure_prefer_prevention_disaster_se <- lmtest::coeftest(exposure_prefer_prevention_disaster_lm, 
                                                                vcov = sandwich::vcovHC, 
                                                                type = "HC1")[,2]

# Put models and robust standard errors together in lists

exposure_lm <- list(exposure_more_relief_covid_lm,
                    exposure_more_relief_disaster_lm,
                    exposure_more_prevention_covid_lm,
                    exposure_more_prevention_disaster_lm,
                    exposure_prefer_prevention_covid_lm,
                    exposure_prefer_prevention_disaster_lm)

exposure_se <- list(exposure_more_relief_covid_se,
                    exposure_more_relief_disaster_se,
                    exposure_more_prevention_covid_se,
                    exposure_more_prevention_disaster_se,
                    exposure_prefer_prevention_covid_se,
                    exposure_prefer_prevention_disaster_se)

# Clean covariate labels

covariate_labels <- c("Age: 40-59", "Age: 60 or older", 
                      "Race: Black", "Race: Asian", "Race: Multiple and Other",
                      "Gender: Female", "Gender: Other", 
                      "Party: Independent", "Party: Democrat", 
                      "Exposure: Covid", "Exposure: Disaster")

# Table A.19

stargazer::stargazer(exposure_lm,
                     se = exposure_se,
                     covariate.labels = covariate_labels,
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     dep.var.labels = c("More on relief (five-point scale, 0-1)","More on prev. (five-point scale, 0-1)","Prioritize prev. (binary)"),
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/exposure-tests.tex",
                     label = "exposure-tests")
```

\input{tables/exposure-tests.tex}

\pagebreak

See Stata file for testing whether the effect of exposure is significantly different across outcomes, and across exposure types within the same outcomes (discussed in section A.7 and A.8 of the online appendix).

The code below generates Table A.20 in the online appendix.

```{r opentext-lm-2, results = 'hide'}

# In this robustness check, we estimate the linear relationship between mentions and preferences but use "only" codings (the respondent only mentions relief, only mentions prevention, or mentions both) and allow all three predictors to enter regression at the same time.

# Regress support for spending more on relief on all three predictors

opentext_lm_more_relief_only <- lm(more_relief ~ mention_only_relief + mention_only_prev + mention_both, data = opentext_responses)
opentext_se_more_relief_only <- sqrt(diag(sandwich::vcovHC(opentext_lm_more_relief_only, type = "HC1")))

# Regress support for spending more on prevention on all three predictors

opentext_lm_more_prevention_only <- lm(more_prevention ~ mention_only_relief + mention_only_prev + mention_both, data = opentext_responses)
opentext_se_more_prevention_only <- sqrt(diag(sandwich::vcovHC(opentext_lm_more_prevention_only, type = "HC1")))

# Combine the regression model objects and vectors containing robust standard errors separately into lists

opentext_lm_2 <- list(opentext_lm_more_relief_only,
                      opentext_lm_more_prevention_only)

opentext_se_2 <- list(opentext_se_more_relief_only,
                      opentext_se_more_prevention_only)

# Generate tables (using stargazer)

# Table A.20

stargazer::stargazer(opentext_lm_2,
                     se = opentext_se_2,
                     dep.var.caption = "",
                     dep.var.labels = c("More on relief", "More on prevention"),
                     covariate.labels = c("Mention only relief", "Mention only prevention", "Mention both"),
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/opentext-validation-only.tex",
                     label = "opentext-validation-only")
```

\input{tables/opentext-validation-only.tex}

\pagebreak

The code below generates Table A.21 in the online appendix.

```{r open-text-lm-1-with-controls, results = 'hide'}

# Run the same regression as for Table 1 but with controls for demographics and party ID

opentext_lm_more_relief_controls <- lm(more_relief ~ age + gender + race + party + mention_relief, data = opentext_responses)
opentext_se_more_relief_controls <- sqrt(diag(sandwich::vcovHC(opentext_lm_more_relief_controls, type = "HC1")))

opentext_lm_more_prevention_controls <- lm(more_prevention ~ age + gender + race + party + mention_prev, data = opentext_responses)
opentext_se_more_prevention_controls <- sqrt(diag(sandwich::vcovHC(opentext_lm_more_prevention_controls, type = "HC1")))

opentext_lm_prefer_prevention_mention_relief_controls <- lm(prefer_prevention ~ age + gender + race + party + mention_relief, data = opentext_responses)
opentext_se_prefer_prevention_mention_relief_controls <- sqrt(diag(sandwich::vcovHC(opentext_lm_prefer_prevention_mention_relief_controls, type = "HC1")))

opentext_lm_prefer_prevention_mention_prevention_controls <- lm(prefer_prevention ~ age + gender + race + party + mention_prev, data = opentext_responses)
opentext_se_prefer_prevention_mention_prevention_controls <- sqrt(diag(sandwich::vcovHC(opentext_lm_prefer_prevention_mention_prevention_controls, type = "HC1")))

# Combine the regression model objects and vectors containing robust standard errors separately into lists

opentext_lm_1_controls <- list(opentext_lm_more_relief_controls,
                               opentext_lm_more_prevention_controls,
                               opentext_lm_prefer_prevention_mention_relief_controls,
                               opentext_lm_prefer_prevention_mention_prevention_controls
                               )

opentext_se_1_controls <- list(opentext_se_more_relief_controls,
                               opentext_se_more_prevention_controls,
                               opentext_se_prefer_prevention_mention_relief_controls,
                               opentext_se_prefer_prevention_mention_prevention_controls
                               )

# Table A.21

stargazer::stargazer(opentext_lm_1_controls,
                     se = opentext_se_1_controls,
                     dep.var.caption = "",
                     dep.var.labels = c("More on relief", "More on prevention", "Prioritize prevention"),
                     covariate.labels = c("Age: 40-59", "Age: 60 or older", 
                                          "Gender: Female", "Gender: Other",
                                          "Race: Black", "Race: Asian", "Race: Multiple and Other",
                                          "Party: Independent", "Party: Democrat", 
                                          "Mention relief", "Mention prevention"),
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/opentext-validation-controls.tex",
                     label = "opentext-validation-controls")
```

\input{tables/opentext-validation-controls.tex}

\pagebreak

The code below generates Table A.22 in the online appendix.

```{r opentext-glm-1, results = 'hide'}

# In this robustness check, we estimate the relationship between mentions and preferences but use logistic regression.
# Note that for some reason, coeftest output with glm models looks different, so to extract standard errors we take the second column only

# Regress support for spending more on relief on whether the respondent mentioned relief in open text

opentext_glm_more_relief <- glm(more_relief ~ mention_relief, data = opentext_responses_binary, binomial(link = "logit"))
opentext_glm_se_more_relief <- lmtest::coeftest(opentext_glm_more_relief, 
                                                  vcov = sandwich::vcovHC, 
                                                  type = "HC1")[,2]

# Run the same regression with controls for demographics and party ID

opentext_glm_more_relief_controls <- glm(more_relief ~ age + gender + race + party + mention_relief, data = opentext_responses_binary, binomial(link = "logit"))
opentext_glm_se_more_relief_controls <- lmtest::coeftest(opentext_glm_more_relief_controls, 
                                                               vcov = sandwich::vcovHC, 
                                                               type = "HC1")[,2]

# Regress support for spending more on prevention on whether the respondent mentioned prevention in open text

opentext_glm_more_prevention <- glm(more_prevention ~ mention_prev, data = opentext_responses_binary, binomial(link = "logit"))
opentext_glm_se_more_prevention <- lmtest::coeftest(opentext_glm_more_prevention, 
                                                          vcov = sandwich::vcovHC, 
                                                          type = "HC1")[,2]

# Run the same regression with controls for demographics and party ID

opentext_glm_more_prevention_controls <- glm(more_prevention ~ age + gender + race + party + mention_prev, data = opentext_responses_binary, binomial(link = "logit"))
opentext_glm_se_more_prevention_controls <- lmtest::coeftest(opentext_glm_more_prevention_controls, 
                                                                   vcov = sandwich::vcovHC, 
                                                                   type = "HC1")[,2]

# Regress support for prioritizing prevention over relief on whether the respondent mentioned relief in open text

opentext_glm_prefer_prevention_mention_relief <- glm(prefer_prevention ~ mention_relief, data = opentext_responses_binary, binomial(link = "logit"))
opentext_glm_se_prefer_prevention_mention_relief <- lmtest::coeftest(opentext_glm_prefer_prevention_mention_relief, 
                                                                           vcov = sandwich::vcovHC, 
                                                                           type = "HC1")[,2]

# Run the same regression with controls for demographics and party ID

opentext_glm_prefer_prevention_mention_relief_controls <- glm(prefer_prevention ~ age + gender + race + party + mention_relief, data = opentext_responses_binary, binomial(link = "logit"))
opentext_glm_se_prefer_prevention_mention_relief_controls <- lmtest::coeftest(opentext_glm_prefer_prevention_mention_relief_controls, 
                                                                                    vcov = sandwich::vcovHC, 
                                                                                    type = "HC1")[,2]

# Regress support for prioritizing prevention over relief on whether the respondent mentioned prevention in open text

opentext_glm_prefer_prevention_mention_prevention <- glm(prefer_prevention ~ mention_prev, data = opentext_responses_binary, binomial(link = "logit"))
opentext_glm_se_prefer_prevention_mention_prevention <- lmtest::coeftest(opentext_glm_prefer_prevention_mention_prevention, 
                                                                               vcov = sandwich::vcovHC, 
                                                                               type = "HC1")[,2]

# Run the same regression with controls for demographics and party ID

opentext_glm_prefer_prevention_mention_prevention_controls <- glm(prefer_prevention ~ age + gender + race + party + mention_prev, data = opentext_responses_binary, binomial(link = "logit"))
opentext_glm_se_prefer_prevention_mention_prevention_controls <- lmtest::coeftest(opentext_glm_prefer_prevention_mention_prevention_controls, 
                                                                                        vcov = sandwich::vcovHC, 
                                                                                        type = "HC1")[,2]

# Combine the regression model objects and vectors containing robust standard errors separately into lists


opentext_lm_3 <- list(opentext_glm_more_relief_controls,
                      opentext_glm_more_relief,
                      opentext_glm_more_prevention_controls,
                      opentext_glm_more_prevention,
                      opentext_glm_prefer_prevention_mention_relief_controls,
                      opentext_glm_prefer_prevention_mention_relief,
                      opentext_glm_prefer_prevention_mention_prevention_controls,
                      opentext_glm_prefer_prevention_mention_prevention)

opentext_se_3 <- list(opentext_glm_se_more_relief_controls,
                      opentext_glm_se_more_relief,
                      opentext_glm_se_more_prevention_controls,
                      opentext_glm_se_more_prevention,
                      opentext_glm_se_prefer_prevention_mention_relief_controls,
                      opentext_glm_se_prefer_prevention_mention_relief,
                      opentext_glm_se_prefer_prevention_mention_prevention_controls,
                      opentext_glm_se_prefer_prevention_mention_prevention)

# Generate table (using stargazer)

# Table A.22

stargazer::stargazer(opentext_lm_3,
                     se = opentext_se_3,
                     dep.var.caption = "",
                     dep.var.labels = c("More on relief", "More on prevention", "Prioritize prevention"),
                     covariate.labels = c("Age: 40-59", "Age: 60 or older", 
                                          "Gender: Female", "Gender: Other",
                                          "Race: Black", "Race: Asian", "Race: Multiple and Other",
                                          "Party: Independent", "Party: Democrat", 
                                          "Mention relief", "Mention prevention"),
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     no.space = TRUE,
                     df = FALSE,
                     header = FALSE,
                     font.size = "tiny",
                     style = "apsr",
                     omit.stat = c("f"),
                     out = "tables/opentext-validation-glm.tex",
                     label = "opentext-validation-glm")
```

\input{tables/opentext-validation-glm.tex}

\pagebreak

## Additional Statistics

The code below produces the correlation coefficients between pandemic-specific and non-pandemic-specific responses) reported in section A.7 of the online appendix (see p. 21-22).

```{r pandemic-correlations}

# Drop missing values and estimate correlation

print("Correlation: More pandemic relief and more relief")
non_missing_relief <- survey_1 %>% select(ResponseId, more_relief, more_pandemic_relief) %>% filter(!is.na(more_pandemic_relief), !is.na(more_relief))
cor(non_missing_relief$more_pandemic_relief, non_missing_relief$more_relief)

print("Correlation: More pandemic prevention and more prevention")
non_missing_prevention <- survey_1 %>% select(ResponseId, more_prevention, more_pandemic_prevention) %>% filter(!is.na(more_pandemic_prevention), !is.na(more_prevention))
cor(non_missing_prevention$more_pandemic_prevention, non_missing_prevention$more_prevention)

print("Correlation: Prioritize pandemic prevention and prioritize prevention")
non_missing_prefer_prevention <- survey_1 %>% select(ResponseId, prefer_prevention, prefer_pandemic_prevention) %>% filter(!is.na(prefer_pandemic_prevention), !is.na(prefer_prevention))
cor(non_missing_prefer_prevention$prefer_pandemic_prevention, non_missing_prefer_prevention$prefer_prevention)
```





